Edit Mesh: select interior faces now detects regions

The previous method only worked in simple cases where faces were
surrounded by non-manifold edges.

Now face regions surrounded by non-manifold edges are marked as interior.
Starting with regions most perpendicular to the surrounding geometry.

Resolves T68401

1 files changed, 313 insertions, 16 deletions

diff --git a/source/blender/editors/mesh/editmesh_select.c b/source/blender/editors/mesh/editmesh_select.c
index d11a223f455..cad9e9a3d06 100644
--- a/source/blender/editors/mesh/editmesh_select.c
+++ b/source/blender/editors/mesh/editmesh_select.c
@@ -31,6 +31,8 @@
 #include "BLI_math_bits.h"
 #include "BLI_rand.h"
 #include "BLI_array.h"
+#include "BLI_heap.h"
+#include "BLI_utildefines_stack.h"
 
 #include "BKE_context.h"
 #include "BKE_report.h"
@@ -2655,39 +2657,334 @@ bool EDBM_mesh_deselect_all_multi(struct bContext *C)
 /* -------------------------------------------------------------------- */
 /** \name Select Interior Faces
  *
- * \note This algorithm is limited to single faces and could be improved, see:
- * https://blender.stackexchange.com/questions/18916
+ * Overview of the algorithm:
+ * - Groups faces surrounded by edges with 3+ faces using them.
+ * - Calculates a cost of each face group comparing it's angle with the faces
+ *   connected to it's non-manifold edges.
+ * - Mark the face group as interior, and mark connected face groups for recalculation.
+ * - Continue to remove the face groups with the highest 'cost'.
+ *
  * \{ */
 
+struct BMFaceLink {
+  struct BMFaceLink *next, *prev;
+  BMFace *face;
+  float area;
+};
+
+static bool bm_interior_loop_filter_fn(const BMLoop *l, void *UNUSED(user_data))
+{
+  if (BM_elem_flag_test(l->e, BM_ELEM_TAG)) {
+    return false;
+  }
+  return true;
+}
+static bool bm_interior_edge_is_manifold_except_face_index(BMEdge *e,
+                                                           int face_index,
+                                                           BMLoop *r_l_pair[2])
+{
+
+  BMLoop *l_iter = e->l;
+  int loop_index = 0;
+  do {
+    BMFace *f = l_iter->f;
+    int i = BM_elem_index_get(f);
+    if (!ELEM(i, -1, face_index)) {
+      if (loop_index == 2) {
+        return false;
+      }
+      r_l_pair[loop_index++] = l_iter;
+    }
+  } while ((l_iter = l_iter->radial_next) != e->l);
+  return (loop_index == 2);
+}
+
+/**
+ * Calculate the cost of the face group.
+ * A higher value means it's more likely to remove first.
+ */
+static float bm_interior_face_group_calc_cost(ListBase *ls, const float *edge_lengths)
+{
+  /* Dividing by the area is important so larger face groups (which will become the outer shell)
+   * aren't detected as having a high cost. */
+  float area = 0.0f;
+  float cost = 0.0f;
+  bool found = false;
+  LISTBASE_FOREACH (struct BMFaceLink *, f_link, ls) {
+    BMFace *f = f_link->face;
+    area += f_link->area;
+    int i = BM_elem_index_get(f);
+    BLI_assert(i != -1);
+    BMLoop *l_iter, *l_first;
+    l_iter = l_first = BM_FACE_FIRST_LOOP(f);
+    do {
+      if (BM_elem_flag_test(l_iter->e, BM_ELEM_TAG)) {
+        float cost_test = 0.0f;
+        int cost_count = 0;
+        /* All other faces. */
+        BMLoop *l_radial_iter = l_iter;
+        do {
+          int i_other = BM_elem_index_get(l_radial_iter->f);
+          if (!ELEM(i_other, -1, i)) {
+            float angle = angle_normalized_v3v3(f->no, l_radial_iter->f->no);
+            /* Ignore face direction since in the case on non-manifold faces connecting edges,
+             * the face flipping may not be meaningful. */
+            if (angle > DEG2RADF(90)) {
+              angle = DEG2RADF(180) - angle;
+            }
+            /* Avoid calculating it inline, pass in pre-calculated edge lengths. */
+#if 0
+            cost_test += BM_edge_calc_length(l_iter->e) * angle;
+#else
+            BLI_assert(edge_lengths[BM_elem_index_get(l_iter->e)] != -1.0f);
+            cost_test += edge_lengths[BM_elem_index_get(l_iter->e)] * angle;
+#endif
+            cost_count += 1;
+          }
+        } while ((l_radial_iter = l_radial_iter->radial_next) != l_iter);
+
+        if (cost_count >= 2) {
+          cost += cost_test;
+          found = true;
+        }
+      }
+    } while ((l_iter = l_iter->next) != l_first);
+  }
+  return found ? cost / area : FLT_MAX;
+}
+
 bool EDBM_select_interior_faces(BMEditMesh *em)
 {
   BMesh *bm = em->bm;
   BMIter iter;
-  BMIter eiter;
-  BMFace *efa;
-  BMEdge *eed;
-  bool ok;
   bool changed = false;
 
-  BM_ITER_MESH (efa, &iter, em->bm, BM_FACES_OF_MESH) {
-    if (BM_elem_flag_test(efa, BM_ELEM_HIDDEN)) {
-      continue;
+  float *edge_lengths = MEM_mallocN(sizeof(*edge_lengths) * bm->totedge, __func__);
+
+  {
+    bool has_nonmanifold = false;
+    BMEdge *e;
+    int i;
+    BM_ITER_MESH_INDEX (e, &iter, bm, BM_EDGES_OF_MESH, i) {
+      const bool is_over = BM_edge_face_count_is_over(e, 2);
+      if (is_over) {
+        BM_elem_flag_enable(e, BM_ELEM_TAG);
+        has_nonmanifold = true;
+        edge_lengths[i] = BM_edge_calc_length(e);
+      }
+      else {
+        BM_elem_flag_disable(e, BM_ELEM_TAG);
+        edge_lengths[i] = -1.0;
+      }
+
+      BM_elem_index_set(e, i); /* set_inline */
     }
+    bm->elem_index_dirty &= ~BM_EDGE;
 
-    ok = true;
-    BM_ITER_ELEM (eed, &eiter, efa, BM_EDGES_OF_FACE) {
-      if (!BM_edge_face_count_is_over(eed, 2)) {
-        ok = false;
+    if (has_nonmanifold == false) {
+      MEM_freeN(edge_lengths);
+      return false;
+    }
+  }
+
+  /* group vars */
+  int *fgroup_array;
+  int(*fgroup_index)[2];
+  int fgroup_len;
+
+  fgroup_array = MEM_mallocN(sizeof(*fgroup_array) * bm->totface, __func__);
+  fgroup_len = BM_mesh_calc_face_groups(
+      bm, fgroup_array, &fgroup_index, bm_interior_loop_filter_fn, NULL, 0, BM_EDGE);
+
+  int *fgroup_recalc_stack = MEM_mallocN(sizeof(*fgroup_recalc_stack) * fgroup_len, __func__);
+  STACK_DECLARE(fgroup_recalc_stack);
+  STACK_INIT(fgroup_recalc_stack, fgroup_len);
+
+  BM_mesh_elem_table_ensure(bm, BM_FACE);
+
+  {
+    BMFace *f;
+    BM_ITER_MESH (f, &iter, em->bm, BM_FACES_OF_MESH) {
+      BM_elem_index_set(f, -1); /* set_dirty! */
+    }
+  }
+  bm->elem_index_dirty |= BM_FACE;
+
+  ListBase *fgroup_listbase = MEM_callocN(sizeof(*fgroup_listbase) * fgroup_len, __func__);
+  struct BMFaceLink *f_link_array = MEM_callocN(sizeof(*f_link_array) * bm->totface, __func__);
+
+  for (int i = 0; i < fgroup_len; i++) {
+    const int fg_sta = fgroup_index[i][0];
+    const int fg_len = fgroup_index[i][1];
+    for (int j = 0; j < fg_len; j++) {
+      const int face_index = fgroup_array[fg_sta + j];
+      BMFace *f = BM_face_at_index(bm, face_index);
+      BM_elem_index_set(f, i);
+
+      struct BMFaceLink *f_link = &f_link_array[face_index];
+      f_link->face = f;
+      f_link->area = BM_face_calc_area(f);
+      BLI_addtail(&fgroup_listbase[i], f_link);
+    }
+  }
+
+  MEM_freeN(fgroup_array);
+  MEM_freeN(fgroup_index);
+
+  Heap *fgroup_heap = BLI_heap_new_ex(fgroup_len);
+  HeapNode **fgroup_table = MEM_mallocN(sizeof(*fgroup_table) * fgroup_len, __func__);
+  bool *fgroup_dirty = MEM_callocN(sizeof(*fgroup_dirty) * fgroup_len, __func__);
+
+  for (int i = 0; i < fgroup_len; i++) {
+    const float cost = bm_interior_face_group_calc_cost(&fgroup_listbase[i], edge_lengths);
+    if (cost != FLT_MAX) {
+      fgroup_table[i] = BLI_heap_insert(fgroup_heap, -cost, POINTER_FROM_INT(i));
+    }
+    else {
+      fgroup_table[i] = NULL;
+    }
+  }
+
+  /* Avoid re-running cost calculations for large face-groups which will end up forming the
+   * outer shell and not be considered interior.
+   * As these face groups become increasingly bigger - their chance of being considered
+   * interior reduces as does the time to calculate their cost.
+   *
+   * This delays recalculating them until they are considered can dates to remove
+   * which becomes less and less likely as they increase in area. */
+
+#define USE_DELAY_FACE_GROUP_COST_CALC
+
+  while (true) {
+
+#if defined(USE_DELAY_FACE_GROUP_COST_CALC)
+    while (!BLI_heap_is_empty(fgroup_heap)) {
+      HeapNode *node_min = BLI_heap_top(fgroup_heap);
+      const int i = POINTER_AS_INT(BLI_heap_node_ptr(node_min));
+      if (fgroup_dirty[i]) {
+        const float cost = bm_interior_face_group_calc_cost(&fgroup_listbase[i], edge_lengths);
+        if (cost != FLT_MAX) {
+          /* The cost may have improves (we may be able to skip this),
+           * however the cost should _never_ make this a choice. */
+          BLI_assert(-BLI_heap_node_value(node_min) >= cost);
+          BLI_heap_node_value_update(fgroup_heap, fgroup_table[i], -cost);
+        }
+        else {
+          BLI_heap_remove(fgroup_heap, fgroup_table[i]);
+          fgroup_table[i] = NULL;
+        }
+        fgroup_dirty[i] = false;
+      }
+      else {
         break;
       }
     }
+#endif
 
-    if (ok) {
-      BM_face_select_set(bm, efa, true);
-      changed = true;
+    if (BLI_heap_is_empty(fgroup_heap)) {
+      break;
     }
+
+    const int i_min = POINTER_AS_INT(BLI_heap_pop_min(fgroup_heap));
+    BLI_assert(fgroup_table[i_min] != NULL);
+    BLI_assert(fgroup_dirty[i_min] == false);
+    fgroup_table[i_min] = NULL;
+    changed = true;
+
+    struct BMFaceLink *f_link;
+    while ((f_link = BLI_pophead(&fgroup_listbase[i_min]))) {
+      BMFace *f = f_link->face;
+      BM_face_select_set(bm, f, true);
+      BM_elem_index_set(f, -1); /* set-dirty */
+
+      BMLoop *l_iter, *l_first;
+
+      /* Loop over edges face edges, merging groups which are no longer separated
+       * by non-manifold edges (when manifold check ignores faces from this group). */
+      l_iter = l_first = BM_FACE_FIRST_LOOP(f);
+      do {
+        BMLoop *l_pair[2];
+        if (bm_interior_edge_is_manifold_except_face_index(l_iter->e, i_min, l_pair)) {
+          BM_elem_flag_disable(l_iter->e, BM_ELEM_TAG);
+
+          int i_a = BM_elem_index_get(l_pair[0]->f);
+          int i_b = BM_elem_index_get(l_pair[1]->f);
+          if (i_a != i_b) {
+            /* Only for predictable results that don't depend on the order of radial loops,
+             * not essential.  */
+            if (i_a > i_b) {
+              SWAP(int, i_a, i_b);
+            }
+
+            /* Merge the the groups. */
+            LISTBASE_FOREACH (LinkData *, n, &fgroup_listbase[i_b]) {
+              BMFace *f_iter = n->data;
+              BM_elem_index_set(f_iter, i_a);
+            }
+            BLI_movelisttolist(&fgroup_listbase[i_a], &fgroup_listbase[i_b]);
+
+            /* This may have been added to 'fgroup_recalc_stack', instead of removing it,
+             * just check the heap node isn't NULL before recalculating. */
+            BLI_heap_remove(fgroup_heap, fgroup_table[i_b]);
+            fgroup_table[i_b] = NULL;
+            /* Keep the dirty flag as-is for 'i_b', because it may be in the 'fgroup_recalc_stack'
+             * and we don't want to add it again.
+             * Instead rely on the 'fgroup_table[i_b]' being NULL as a secondary check. */
+
+            if (fgroup_dirty[i_a] == false) {
+              BLI_assert(fgroup_table[i_a] != NULL);
+              STACK_PUSH(fgroup_recalc_stack, i_a);
+              fgroup_dirty[i_a] = true;
+            }
+          }
+        }
+
+        /* Mark all connected groups for re-calculation. */
+        BMLoop *l_radial_iter = l_iter->radial_next;
+        if (l_radial_iter != l_iter) {
+          do {
+            int i_other = BM_elem_index_get(l_radial_iter->f);
+            if (!ELEM(i_other, -1, i_min)) {
+              if ((fgroup_table[i_other] != NULL) && (fgroup_dirty[i_other] == false)) {
+#if !defined(USE_DELAY_FACE_GROUP_COST_CALC)
+                STACK_PUSH(fgroup_recalc_stack, i_other);
+#endif
+                fgroup_dirty[i_other] = true;
+              }
+            }
+          } while ((l_radial_iter = l_radial_iter->radial_next) != l_iter);
+        }
+
+      } while ((l_iter = l_iter->next) != l_first);
+    }
+
+    for (int index = 0; index < STACK_SIZE(fgroup_recalc_stack); index++) {
+      const int i = fgroup_recalc_stack[index];
+      if (fgroup_table[i] != NULL && fgroup_dirty[i] == true) {
+        /* First update edge tags. */
+        const float cost = bm_interior_face_group_calc_cost(&fgroup_listbase[i], edge_lengths);
+        if (cost != FLT_MAX) {
+          BLI_heap_node_value_update(fgroup_heap, fgroup_table[i], -cost);
+        }
+        else {
+          BLI_heap_remove(fgroup_heap, fgroup_table[i]);
+          fgroup_table[i] = NULL;
+        }
+      }
+      fgroup_dirty[i] = false;
+    }
+    STACK_CLEAR(fgroup_recalc_stack);
   }
 
+  MEM_freeN(edge_lengths);
+  MEM_freeN(f_link_array);
+  MEM_freeN(fgroup_listbase);
+  MEM_freeN(fgroup_recalc_stack);
+  MEM_freeN(fgroup_table);
+  MEM_freeN(fgroup_dirty);
+
+  BLI_heap_free(fgroup_heap, NULL);
+
   return changed;
 }