Shrinkwrap: improve triangle boundary stability in Target Normal Project.

Rewrite the checks for determining if the solution is actually within
the triangle to fix stability issues when the correct solution is on
an edge, and step is very small, i.e. the solution is already very
close. Also, comment more clearly what is happening geometrically.

This should fix problems when vertices that should project exactly
onto an edge actually miss, resulting in weird spikes. This made
Target Normal Project unusable for the voxel remesher.

1 files changed, 39 insertions, 34 deletions

diff --git a/source/blender/blenkernel/intern/shrinkwrap.c b/source/blender/blenkernel/intern/shrinkwrap.c
index c6cac2057d6..a2699e767e9 100644
--- a/source/blender/blenkernel/intern/shrinkwrap.c
+++ b/source/blender/blenkernel/intern/shrinkwrap.c
@@ -794,54 +794,59 @@ static bool target_project_tri_correct(void *UNUSED(userdata),
                                        float x_next[3])
 {
   /* Insignificant correction threshold */
-  const float epsilon = 1e-6f;
-  const float dir_epsilon = 0.05f;
+  const float epsilon = 1e-5f;
+  /* Dot product threshold for checking if step is 'clearly' pointing outside. */
+  const float dir_epsilon = 0.5f;
   bool fixed = false, locked = false;
 
-  /* Weight 0 and 1 boundary check. */
-  for (int i = 0; i < 2; i++) {
-    if (step[i] > x[i]) {
-      if (step[i] > dir_epsilon * fabsf(step[1 - i])) {
-        /* Abort if the solution is clearly outside the domain. */
-        if (x[i] < epsilon) {
-          return false;
-        }
+  /* The barycentric coordinate domain is a triangle bounded by
+   * the X and Y axes, plus the x+y=1 diagonal. First, clamp the
+   * movement against the diagonal. Note that step is subtracted. */
+  float sum = x[0] + x[1];
+  float sstep = -(step[0] + step[1]);
 
-        /* Scale a significant step down to arrive at the boundary. */
-        mul_v3_fl(step, x[i] / step[i]);
-        fixed = true;
-      }
-      else {
-        /* Reset precision errors to stay at the boundary. */
-        step[i] = x[i];
-        fixed = locked = true;
-      }
-    }
-  }
+  if (sum + sstep > 1.0f) {
+    float ldist = 1.0f - sum;
 
-  /* Weight 2 boundary check. */
-  float sum = x[0] + x[1];
-  float sstep = step[0] + step[1];
+    /* If already at the boundary, slide along it. */
+    if (ldist < epsilon * (float)M_SQRT2) {
+      float step_len = len_v2(step);
 
-  if (sum - sstep > 1.0f) {
-    if (sstep < -dir_epsilon * (fabsf(step[0]) + fabsf(step[1]))) {
       /* Abort if the solution is clearly outside the domain. */
-      if (sum > 1.0f - epsilon) {
+      if (step_len > epsilon && sstep > step_len * dir_epsilon * (float)M_SQRT2) {
         return false;
       }
 
+      /* Project the new position onto the diagonal. */
+      add_v2_fl(step, (sum + sstep - 1.0f) * 0.5f);
+      fixed = locked = true;
+    }
+    else {
       /* Scale a significant step down to arrive at the boundary. */
-      mul_v3_fl(step, (1.0f - sum) / -sstep);
+      mul_v3_fl(step, ldist / sstep);
       fixed = true;
     }
-    else {
-      /* Reset precision errors to stay at the boundary. */
-      if (locked) {
-        step[0] = step[1] = 0.0f;
+  }
+
+  /* Weight 0 and 1 boundary checks - along axis. */
+  for (int i = 0; i < 2; i++) {
+    if (step[i] > x[i]) {
+      /* If already at the boundary, slide along it. */
+      if (x[i] < epsilon) {
+        float step_len = len_v2(step);
+
+        /* Abort if the solution is clearly outside the domain. */
+        if (step_len > epsilon && (locked || step[i] > step_len * dir_epsilon)) {
+          return false;
+        }
+
+        /* Reset precision errors to stay at the boundary. */
+        step[i] = x[i];
+        fixed = true;
       }
       else {
-        step[0] -= 0.5f * sstep;
-        step[1] = -step[0];
+        /* Scale a significant step down to arrive at the boundary. */
+        mul_v3_fl(step, x[i] / step[i]);
         fixed = true;
       }
     }