BMesh: add utility to calculate normal from a vertex cloud

Extract from BM_verts_sort_radial_plane & simplify.

source/blender/bmesh/intern/bmesh_construct.c

@@ -390,79 +390,13 @@ void BM_verts_sort_radial_plane(BMVert **vert_arr, int len)
   struct SortIntByFloat *vang = BLI_array_alloca(vang, len);
   BMVert **vert_arr_map = BLI_array_alloca(vert_arr_map, len);
 
-  float totv_inv = 1.0f / (float)len;
-  int i = 0;
+  float nor[3], cent[3];
+  int index_tangent = 0;
+  BM_verts_calc_normal_from_cloud_ex(vert_arr, len, nor, cent, &index_tangent);
+  const float *far = vert_arr[index_tangent]->co;
 
-  float cent[3], nor[3];
-
-  const float *far = NULL, *far_cross = NULL;
-
-  float far_vec[3];
-  float far_cross_vec[3];
-  float sign_vec[3]; /* work out if we are pos/neg angle */
-
-  float far_dist_sq, far_dist_max_sq;
-  float far_cross_dist, far_cross_best = 0.0f;
-
-  /* get the center point and collect vector array since we loop over these a lot */
-  zero_v3(cent);
-  for (i = 0; i < len; i++) {
-    madd_v3_v3fl(cent, vert_arr[i]->co, totv_inv);
-  }
-
-  /* find the far point from cent */
-  far_dist_max_sq = 0.0f;
-  for (i = 0; i < len; i++) {
-    far_dist_sq = len_squared_v3v3(vert_arr[i]->co, cent);
-    if (far_dist_sq > far_dist_max_sq || far == NULL) {
-      far = vert_arr[i]->co;
-      far_dist_max_sq = far_dist_sq;
-    }
-  }
-
-  sub_v3_v3v3(far_vec, far, cent);
-  // far_dist = len_v3(far_vec); /* real dist */ /* UNUSED */
-
-  /* --- */
-
-  /* find a point 90deg about to compare with */
-  far_cross_best = 0.0f;
-  for (i = 0; i < len; i++) {
-
-    if (far == vert_arr[i]->co) {
-      continue;
-    }
-
-    sub_v3_v3v3(far_cross_vec, vert_arr[i]->co, cent);
-    far_cross_dist = normalize_v3(far_cross_vec);
-
-    /* more of a weight then a distance */
-    far_cross_dist = (
-        /* First we want to have a value close to zero mapped to 1. */
-        1.0f - fabsf(dot_v3v3(far_vec, far_cross_vec)) *
-                   /* Second  we multiply by the distance
-                    * so points close to the center are not preferred. */
-                   far_cross_dist);
-
-    if (far_cross_dist > far_cross_best || far_cross == NULL) {
-      far_cross = vert_arr[i]->co;
-      far_cross_best = far_cross_dist;
-    }
-  }
-
-  sub_v3_v3v3(far_cross_vec, far_cross, cent);
-
-  /* --- */
-
-  /* now we have 2 vectors we can have a cross product */
-  cross_v3_v3v3(nor, far_vec, far_cross_vec);
-  normalize_v3(nor);
-  cross_v3_v3v3(sign_vec, far_vec, nor); /* this vector should match 'far_cross_vec' closely */
-
-  /* --- */
-
-  /* now calculate every points angle around the normal (signed) */
-  for (i = 0; i < len; i++) {
+  /* Now calculate every points angle around the normal (signed). */
+  for (int i = 0; i < len; i++) {
     vang[i].sort_value = angle_signed_on_axis_v3v3v3_v3(far, cent, vert_arr[i]->co, nor);
     vang[i].data = i;
     vert_arr_map[i] = vert_arr[i];
@@ -473,7 +407,7 @@ void BM_verts_sort_radial_plane(BMVert **vert_arr, int len)
 
   /* --- */
 
-  for (i = 0; i < len; i++) {
+  for (int i = 0; i < len; i++) {
     vert_arr[i] = vert_arr_map[vang[i].data];
   }
 }

source/blender/bmesh/intern/bmesh_polygon.c

@@ -864,6 +864,113 @@ float BM_face_calc_normal_vcos(const BMesh *bm,
   }
 }
 
+/**
+ * Calculate a normal from a vertex cloud.
+ *
+ * \note We could make a higher quality version that takes all vertices into account.
+ * Currently it finds 4 outer most points returning it's normal.
+ */
+void BM_verts_calc_normal_from_cloud_ex(
+    BMVert **varr, int varr_len, float r_normal[3], float r_center[3], int *r_index_tangent)
+{
+  const float varr_len_inv = 1.0f / (float)varr_len;
+
+  /* Get the center point and collect vector array since we loop over these a lot. */
+  float center[3] = {0.0f, 0.0f, 0.0f};
+  for (int i = 0; i < varr_len; i++) {
+    madd_v3_v3fl(center, varr[i]->co, varr_len_inv);
+  }
+
+  /* Find the 'co_a' point from center. */
+  int co_a_index = 0;
+  const float *co_a = NULL;
+  {
+    float dist_sq_max = -1.0f;
+    for (int i = 0; i < varr_len; i++) {
+      const float dist_sq_test = len_squared_v3v3(varr[i]->co, center);
+      if (!(dist_sq_test <= dist_sq_max)) {
+        co_a = varr[i]->co;
+        co_a_index = i;
+        dist_sq_max = dist_sq_test;
+      }
+    }
+  }
+
+  float dir_a[3];
+  sub_v3_v3v3(dir_a, co_a, center);
+  normalize_v3(dir_a);
+
+  const float *co_b = NULL;
+  float dir_b[3] = {0.0f, 0.0f, 0.0f};
+  {
+    float dist_sq_max = -1.0f;
+    for (int i = 0; i < varr_len; i++) {
+      if (varr[i]->co == co_a) {
+        continue;
+      }
+      float dir_test[3];
+      sub_v3_v3v3(dir_test, varr[i]->co, center);
+      project_plane_normalized_v3_v3v3(dir_test, dir_test, dir_a);
+      const float dist_sq_test = len_squared_v3(dir_test);
+      if (!(dist_sq_test <= dist_sq_max)) {
+        co_b = varr[i]->co;
+        dist_sq_max = dist_sq_test;
+        copy_v3_v3(dir_b, dir_test);
+      }
+    }
+  }
+
+  if (varr_len <= 3) {
+    normal_tri_v3(r_normal, center, co_a, co_b);
+    goto finally;
+  }
+
+  normalize_v3(dir_b);
+
+  const float *co_a_opposite = NULL;
+  const float *co_b_opposite = NULL;
+
+  {
+    float dot_a_min = FLT_MAX;
+    float dot_b_min = FLT_MAX;
+    for (int i = 0; i < varr_len; i++) {
+      const float *co_test = varr[i]->co;
+      float dot_test;
+
+      if (co_test != co_a) {
+        dot_test = dot_v3v3(dir_a, co_test);
+        if (dot_test < dot_a_min) {
+          dot_a_min = dot_test;
+          co_a_opposite = co_test;
+        }
+      }
+
+      if (co_test != co_b) {
+        dot_test = dot_v3v3(dir_b, co_test);
+        if (dot_test < dot_b_min) {
+          dot_b_min = dot_test;
+          co_b_opposite = co_test;
+        }
+      }
+    }
+  }
+
+  normal_quad_v3(r_normal, co_a, co_b, co_a_opposite, co_b_opposite);
+
+finally:
+  if (r_center != NULL) {
+    copy_v3_v3(r_center, center);
+  }
+  if (r_index_tangent != NULL) {
+    *r_index_tangent = co_a_index;
+  }
+}
+
+void BM_verts_calc_normal_from_cloud(BMVert **varr, int varr_len, float r_normal[3])
+{
+  BM_verts_calc_normal_from_cloud_ex(varr, varr_len, r_normal, NULL, NULL);
+}
+
 /**
  * Calculates the face subset normal.
  */

source/blender/bmesh/intern/bmesh_polygon.h

@@ -38,6 +38,11 @@ float BM_face_calc_normal_vcos(const BMesh *bm,
                                const BMFace *f,
                                float r_no[3],
                                float const (*vertexCos)[3]) ATTR_NONNULL();
+
+void BM_verts_calc_normal_from_cloud_ex(
+    BMVert **varr, int varr_len, float r_normal[3], float r_center[3], int *r_index_tangent);
+void BM_verts_calc_normal_from_cloud(BMVert **varr, int varr_len, float r_normal[3]);
+
 float BM_face_calc_normal_subset(const BMLoop *l_first, const BMLoop *l_last, float r_no[3])
     ATTR_NONNULL();
 float BM_face_calc_area(const BMFace *f) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();