BKE: parallelize BKE_mesh_calc_edges

`BKE_mesh_calc_edges` was the main performance bottleneck in D9141.
While openvdb only needed ~115ms, calculating the edges afterwards
took ~960ms. Now with some parallelization this is reduced to ~210ms.

Parallelizing `BKE_mesh_calc_edges` is not entirely trivial, because it
has to perform deduplication and some other things that have to happen
in a certain order. Even though the multithreading improves performance
with more threads, there are diminishing returns when too many threads
are used in this function.

The speedup is mainly achieved by having multiple hash tables that are
filled in parallel. The distribution of the edges to hash tables is based on
a hash (that is different from the hash used in the actual hash tables).

I moved the function to C++, because that made it easier for me to
optimize it. Furthermore, I added `BLI_task.hh` which contains some
light tbb wrappers for parallelization.

Reviewers: campbellbarton

Differential Revision: https://developer.blender.org/D9151

5 files changed, 332 insertions, 106 deletions

diff --git a/source/blender/blenkernel/BKE_mesh.h b/source/blender/blenkernel/BKE_mesh.h
index a61e453ec52..fdea26ce730 100644
--- a/source/blender/blenkernel/BKE_mesh.h
+++ b/source/blender/blenkernel/BKE_mesh.h
@@ -673,7 +673,7 @@ void BKE_mesh_strip_loose_edges(struct Mesh *me);
 
 void BKE_mesh_calc_edges_legacy(struct Mesh *me, const bool use_old);
 void BKE_mesh_calc_edges_loose(struct Mesh *mesh);
-void BKE_mesh_calc_edges(struct Mesh *mesh, bool update, const bool select);
+void BKE_mesh_calc_edges(struct Mesh *mesh, bool keep_existing_edges, const bool select_new_edges);
 void BKE_mesh_calc_edges_tessface(struct Mesh *mesh);
 
 /* In DerivedMesh.c */
diff --git a/source/blender/blenkernel/CMakeLists.txt b/source/blender/blenkernel/CMakeLists.txt
index 63e98b94852..0fbc8c4c229 100644
--- a/source/blender/blenkernel/CMakeLists.txt
+++ b/source/blender/blenkernel/CMakeLists.txt
@@ -175,6 +175,7 @@ set(SRC
   intern/mesh_runtime.c
   intern/mesh_tangent.c
   intern/mesh_validate.c
+  intern/mesh_validate.cc
   intern/mesh_wrapper.c
   intern/modifier.c
   intern/movieclip.c
@@ -690,6 +691,18 @@ if(WITH_XR_OPENXR)
   add_definitions(-DWITH_XR_OPENXR)
 endif()
 
+if(WITH_TBB)
+  add_definitions(-DWITH_TBB)
+
+  list(APPEND INC_SYS
+    ${TBB_INCLUDE_DIRS}
+  )
+
+  list(APPEND LIB
+    ${TBB_LIBRARIES}
+  )
+endif()
+
 # # Warnings as errors, this is too strict!
 # if(MSVC)
 #    set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} /WX")
diff --git a/source/blender/blenkernel/intern/mesh_validate.c b/source/blender/blenkernel/intern/mesh_validate.c
index af09ef6dae0..5f11476574e 100644
--- a/source/blender/blenkernel/intern/mesh_validate.c
+++ b/source/blender/blenkernel/intern/mesh_validate.c
@@ -1526,111 +1526,6 @@ void BKE_mesh_calc_edges_legacy(Mesh *me, const bool use_old)
   BKE_mesh_strip_loose_faces(me);
 }
 
-/**
- * Calculate edges from polygons
- *
- * \param mesh: The mesh to add edges into
- * \param update: When true create new edges co-exist
- */
-void BKE_mesh_calc_edges(Mesh *mesh, bool update, const bool select)
-{
-  MEdge *med;
-  int i, totpoly = mesh->totpoly;
-  /* Select for newly created meshes which are selected T25595. */
-  const short ed_flag = (ME_EDGEDRAW | ME_EDGERENDER) | (select ? SELECT : 0);
-
-  if (mesh->totedge == 0) {
-    update = false;
-  }
-
-  const unsigned int eh_reserve = max_ii(update ? mesh->totedge : 0,
-                                         BLI_EDGEHASH_SIZE_GUESS_FROM_POLYS(totpoly));
-  EdgeHash *eh = BLI_edgehash_new_ex(__func__, eh_reserve);
-
-  if (update) {
-    /* assume existing edges are valid
-     * useful when adding more faces and generating edges from them */
-    med = mesh->medge;
-    for (i = 0; i < mesh->totedge; i++, med++) {
-      BLI_edgehash_insert(eh, med->v1, med->v2, med);
-    }
-  }
-
-  /* mesh loops (bmesh only) */
-  MPoly *mp;
-  for (mp = mesh->mpoly, i = 0; i < totpoly; mp++, i++) {
-    MLoop *l = &mesh->mloop[mp->loopstart];
-    int j, v_prev = (l + (mp->totloop - 1))->v;
-    for (j = 0; j < mp->totloop; j++, l++) {
-      if (v_prev != l->v) {
-        void **val_p;
-        if (!BLI_edgehash_ensure_p(eh, v_prev, l->v, &val_p)) {
-          *val_p = NULL;
-        }
-      }
-      v_prev = l->v;
-    }
-  }
-
-  const int totedge = BLI_edgehash_len(eh);
-
-  /* write new edges into a temporary CustomData */
-  CustomData edata;
-  CustomData_reset(&edata);
-  CustomData_add_layer(&edata, CD_MEDGE, CD_CALLOC, NULL, totedge);
-
-  med = CustomData_get_layer(&edata, CD_MEDGE);
-  EdgeHashIterator *ehi;
-  for (ehi = BLI_edgehashIterator_new(eh), i = 0; BLI_edgehashIterator_isDone(ehi) == false;
-       BLI_edgehashIterator_step(ehi), ++i, ++med) {
-    MEdge *med_orig;
-    if (update && (med_orig = BLI_edgehashIterator_getValue(ehi))) {
-      *med = *med_orig; /* copy from the original */
-    }
-    else {
-      BLI_edgehashIterator_getKey(ehi, &med->v1, &med->v2);
-      med->flag = ed_flag;
-    }
-
-    /* store the new edge index in the hash value */
-    BLI_edgehashIterator_setValue(ehi, POINTER_FROM_INT(i));
-  }
-  BLI_edgehashIterator_free(ehi);
-
-  if (mesh->totpoly) {
-    /* second pass, iterate through all loops again and assign
-     * the newly created edges to them. */
-    for (mp = mesh->mpoly, i = 0; i < mesh->totpoly; mp++, i++) {
-      MLoop *l = &mesh->mloop[mp->loopstart];
-      MLoop *l_prev = (l + (mp->totloop - 1));
-      int j;
-      for (j = 0; j < mp->totloop; j++, l++) {
-        /* Lookup hashed edge index, if it's valid. */
-        int med_index;
-        if (l_prev->v != l->v) {
-          med_index = POINTER_AS_INT(BLI_edgehash_lookup(eh, l_prev->v, l->v));
-        }
-        else {
-          /* This is an invalid edge; normally this does not happen in Blender, but it can be part
-           * of an imported mesh with invalid geometry. See T76514. */
-          med_index = 0;
-        }
-        l_prev->e = med_index;
-        l_prev = l;
-      }
-    }
-  }
-
-  /* free old CustomData and assign new one */
-  CustomData_free(&mesh->edata, mesh->totedge);
-  mesh->edata = edata;
-  mesh->totedge = totedge;
-
-  mesh->medge = CustomData_get_layer(&mesh->edata, CD_MEDGE);
-
-  BLI_edgehash_free(eh, NULL);
-}
-
 void BKE_mesh_calc_edges_loose(Mesh *mesh)
 {
   MEdge *med = mesh->medge;
diff --git a/source/blender/blenkernel/intern/mesh_validate.cc b/source/blender/blenkernel/intern/mesh_validate.cc
new file mode 100644
index 00000000000..64fddc81124
--- /dev/null
+++ b/source/blender/blenkernel/intern/mesh_validate.cc
@@ -0,0 +1,255 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/** \file
+ * \ingroup bke
+ */
+
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_object_types.h"
+
+#include "BLI_edgehash.h"
+#include "BLI_map.hh"
+#include "BLI_math_base.h"
+#include "BLI_task.hh"
+#include "BLI_threads.h"
+#include "BLI_timeit.hh"
+
+#include "BKE_customdata.h"
+#include "BKE_mesh.h"
+
+namespace blender::bke::calc_edges {
+
+/** This is used to uniquely identify edges in a hash map. */
+struct OrderedEdge {
+  int v_low, v_high;
+
+  OrderedEdge(const int v1, const int v2)
+  {
+    if (v1 < v2) {
+      v_low = v1;
+      v_high = v2;
+    }
+    else {
+      v_low = v2;
+      v_high = v1;
+    }
+  }
+
+  OrderedEdge(const uint v1, const uint v2)
+      : OrderedEdge(static_cast<int>(v1), static_cast<int>(v2))
+  {
+  }
+
+  uint64_t hash() const
+  {
+    return (this->v_low << 8) ^ this->v_high;
+  }
+
+  /** Return a hash value that is likely to be different in the low bits from the normal `hash()`
+   * function. This is necessary to avoid collisions in #BKE_mesh_calc_edges. */
+  uint64_t hash2() const
+  {
+    return this->v_low;
+  }
+
+  friend bool operator==(const OrderedEdge &e1, const OrderedEdge &e2)
+  {
+    BLI_assert(e1.v_low < e1.v_high);
+    BLI_assert(e2.v_low < e2.v_high);
+    return e1.v_low == e2.v_low && e1.v_high == e2.v_high;
+  }
+};
+
+/* The map first contains an edge pointer and later an index. */
+typedef union OrigEdgeOrIndex {
+  const MEdge *original_edge;
+  int index;
+} OrigEdgeOrIndex;
+using EdgeMap = Map<OrderedEdge, OrigEdgeOrIndex>;
+
+static void add_existing_edges_to_hash_maps(Mesh *mesh,
+                                            MutableSpan<EdgeMap> edge_maps,
+                                            uint32_t parallel_mask)
+{
+  /* Assume existing edges are valid. */
+  parallel_for_each(edge_maps, [&](EdgeMap &edge_map) {
+    const int task_index = &edge_map - &edge_maps[0];
+    for (const MEdge &edge : Span(mesh->medge, mesh->totedge)) {
+      OrderedEdge ordered_edge{edge.v1, edge.v2};
+      /* Only add the edge when it belongs into this map. */
+      if (task_index == (parallel_mask & ordered_edge.hash2())) {
+        edge_map.add_new(ordered_edge, {&edge});
+      }
+    }
+  });
+}
+
+static void add_polygon_edges_to_hash_maps(Mesh *mesh,
+                                           MutableSpan<EdgeMap> edge_maps,
+                                           uint32_t parallel_mask)
+{
+  const Span<MLoop> loops{mesh->mloop, mesh->totloop};
+  parallel_for_each(edge_maps, [&](EdgeMap &edge_map) {
+    const int task_index = &edge_map - &edge_maps[0];
+    for (const MPoly &poly : Span(mesh->mpoly, mesh->totpoly)) {
+      Span<MLoop> poly_loops = loops.slice(poly.loopstart, poly.totloop);
+      const MLoop *prev_loop = &poly_loops.last();
+      for (const MLoop &next_loop : poly_loops) {
+        /* Can only be the same when the mesh data is invalid. */
+        if (prev_loop->v != next_loop.v) {
+          OrderedEdge ordered_edge{prev_loop->v, next_loop.v};
+          /* Only add the edge when it belongs into this map. */
+          if (task_index == (parallel_mask & ordered_edge.hash2())) {
+            edge_map.lookup_or_add(ordered_edge, {nullptr});
+          }
+        }
+        prev_loop = &next_loop;
+      }
+    }
+  });
+}
+
+static void serialize_and_initialize_deduplicated_edges(MutableSpan<EdgeMap> edge_maps,
+                                                        MutableSpan<MEdge> new_edges,
+                                                        short new_edge_flag)
+{
+  /* All edges are distributed in the hash tables now. They have to be serialized into a single
+   * array below. To be able to parallelize this, we have to compute edge index offsets for each
+   * map. */
+  Array<int> edge_index_offsets(edge_maps.size());
+  edge_index_offsets[0] = 0;
+  for (const int i : IndexRange(edge_maps.size() - 1)) {
+    edge_index_offsets[i + 1] = edge_index_offsets[i] + edge_maps[i].size();
+  }
+
+  parallel_for_each(edge_maps, [&](EdgeMap &edge_map) {
+    const int task_index = &edge_map - &edge_maps[0];
+
+    int new_edge_index = edge_index_offsets[task_index];
+    for (EdgeMap::MutableItem item : edge_map.items()) {
+      MEdge &new_edge = new_edges[new_edge_index];
+      const MEdge *orig_edge = item.value.original_edge;
+      if (orig_edge != nullptr) {
+        /* Copy values from original edge. */
+        new_edge = *orig_edge;
+      }
+      else {
+        /* Initialize new edge. */
+        new_edge.v1 = item.key.v_low;
+        new_edge.v2 = item.key.v_high;
+        new_edge.flag = new_edge_flag;
+      }
+      item.value.index = new_edge_index;
+      new_edge_index++;
+    }
+  });
+}
+
+static void update_edge_indices_in_poly_loops(Mesh *mesh,
+                                              Span<EdgeMap> edge_maps,
+                                              uint32_t parallel_mask)
+{
+  const MutableSpan<MLoop> loops{mesh->mloop, mesh->totloop};
+  parallel_for(IndexRange(mesh->totpoly), 100, [&](IndexRange range) {
+    for (const int poly_index : range) {
+      MPoly &poly = mesh->mpoly[poly_index];
+      MutableSpan<MLoop> poly_loops = loops.slice(poly.loopstart, poly.totloop);
+
+      MLoop *prev_loop = &poly_loops.last();
+      for (MLoop &next_loop : poly_loops) {
+        int edge_index;
+        if (prev_loop->v != next_loop.v) {
+          OrderedEdge ordered_edge{prev_loop->v, next_loop.v};
+          /* Double lookup: First find the map that contains the edge, then lookup the edge. */
+          const EdgeMap &edge_map = edge_maps[parallel_mask & ordered_edge.hash2()];
+          edge_index = edge_map.lookup(ordered_edge).index;
+        }
+        else {
+          /* This is an invalid edge; normally this does not happen in Blender,
+           * but it can be part of an imported mesh with invalid geometry. See
+           * T76514. */
+          edge_index = 0;
+        }
+        prev_loop->e = edge_index;
+        prev_loop = &next_loop;
+      }
+    }
+  });
+}
+
+static int get_parallel_maps_count(const Mesh *mesh)
+{
+  /* Don't use parallelization when the mesh is small. */
+  if (mesh->totpoly < 1000) {
+    return 1;
+  }
+  /* Use at most 8 separate hash tables. Using more threads has diminishing returns. These threads
+   * can better do something more useful instead. */
+  const int system_thread_count = BLI_system_thread_count();
+  return power_of_2_min_i(std::min(8, system_thread_count));
+}
+
+}  // namespace blender::bke::calc_edges
+
+/**
+ * Calculate edges from polygons.
+ */
+void BKE_mesh_calc_edges(Mesh *mesh, bool keep_existing_edges, const bool select_new_edges)
+{
+  using namespace blender;
+  using namespace blender::bke;
+  using namespace blender::bke::calc_edges;
+
+  /* Parallelization is achieved by having multiple hash tables for different subsets of edges.
+   * Each edge is assigned to one of the hash maps based on the lower bits of a hash value. */
+  const int parallel_maps = get_parallel_maps_count(mesh);
+  BLI_assert(is_power_of_2_i(parallel_maps));
+  const uint32_t parallel_mask = static_cast<uint32_t>(parallel_maps) - 1;
+  Array<EdgeMap> edge_maps(parallel_maps);
+
+  /* Reserve memory in hash tables. */
+  const int totedge_guess = std::max(keep_existing_edges ? mesh->totedge : 0, mesh->totpoly * 2);
+  parallel_for_each(edge_maps,
+                    [&](EdgeMap &edge_map) { edge_map.reserve(totedge_guess / parallel_maps); });
+
+  /* Add all edges. */
+  if (keep_existing_edges) {
+    calc_edges::add_existing_edges_to_hash_maps(mesh, edge_maps, parallel_mask);
+  }
+  calc_edges::add_polygon_edges_to_hash_maps(mesh, edge_maps, parallel_mask);
+
+  /* Compute total number of edges. */
+  int new_totedge = 0;
+  for (EdgeMap &edge_map : edge_maps) {
+    new_totedge += edge_map.size();
+  }
+
+  /* Create new edges. */
+  MutableSpan<MEdge> new_edges{
+      static_cast<MEdge *>(MEM_calloc_arrayN(new_totedge, sizeof(MEdge), __func__)), new_totedge};
+  const short new_edge_flag = (ME_EDGEDRAW | ME_EDGERENDER) | (select_new_edges ? SELECT : 0);
+  calc_edges::serialize_and_initialize_deduplicated_edges(edge_maps, new_edges, new_edge_flag);
+  calc_edges::update_edge_indices_in_poly_loops(mesh, edge_maps, parallel_mask);
+
+  /* Free old CustomData and assign new one. */
+  CustomData_free(&mesh->edata, mesh->totedge);
+  CustomData_reset(&mesh->edata);
+  CustomData_add_layer(&mesh->edata, CD_MEDGE, CD_ASSIGN, new_edges.data(), new_totedge);
+  mesh->totedge = new_totedge;
+  mesh->medge = new_edges.data();
+}
diff --git a/source/blender/blenlib/BLI_task.hh b/source/blender/blenlib/BLI_task.hh
new file mode 100644
index 00000000000..7e60f271e9d
--- /dev/null
+++ b/source/blender/blenlib/BLI_task.hh
@@ -0,0 +1,63 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#ifdef WITH_TBB
+/* Quiet top level deprecation message, unrelated to API usage here. */
+#  define TBB_SUPPRESS_DEPRECATED_MESSAGES 1
+#  include <tbb/tbb.h>
+#endif
+
+#include "BLI_index_range.hh"
+#include "BLI_utildefines.h"
+
+namespace blender {
+
+template<typename Range, typename Function>
+void parallel_for_each(Range &range, const Function &function)
+{
+#ifdef WITH_TBB
+  tbb::parallel_for_each(range, function);
+#else
+  for (auto &value : range) {
+    function(value);
+  }
+#endif
+}
+
+template<typename Function>
+void parallel_for(IndexRange range, int64_t grain_size, const Function &function)
+{
+  if (range.size() == 0) {
+    return;
+  }
+#ifdef WITH_TBB
+  tbb::parallel_for(tbb::blocked_range<int64_t>(range.first(), range.one_after_last(), grain_size),
+                    [&](const tbb::blocked_range<int64_t> &subrange) {
+                      function(IndexRange(subrange.begin(), subrange.size()));
+                    });
+#else
+  UNUSED_VARS(grain_size);
+  function(range);
+#endif
+}
+
+}  // namespace blender