Shrinkwrap: improve numerical stability of Target Normal Project.

* Add proper adjustment for scale in the solver epsilon computation.
* Run at least one full iteration of the solver, even if the initial
  state meets the epsilon requirement.
* When applying offset, blend normal into the offset direction
  as the initial point moves very close to the target mesh.

Also random improvements to debug trace output in the console.

1 files changed, 3 insertions, 3 deletions

diff --git a/source/blender/blenlib/intern/math_solvers.c b/source/blender/blenlib/intern/math_solvers.c
index a1c3d16a404..235589abdab 100644
--- a/source/blender/blenlib/intern/math_solvers.c
+++ b/source/blender/blenlib/intern/math_solvers.c
@@ -215,10 +215,10 @@ bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta,
   fdeltav = len_squared_v3(fdelta);
 
   if (trace) {
-    printf("START (%g, %g, %g) %g\n", x[0], x[1], x[2], fdeltav);
+    printf("START (%g, %g, %g) %g %g\n", x[0], x[1], x[2], fdeltav, epsilon);
   }
 
-  for (int i = 0; i < max_iterations && fdeltav > epsilon; i++) {
+  for (int i = 0; i == 0 || (i < max_iterations && fdeltav > epsilon); i++) {
     /* Newton's method step. */
     func_jacobian(userdata, x, jacobian);
 
@@ -248,7 +248,7 @@ bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta,
     }
 
     /* Line search correction. */
-    while (next_fdeltav > fdeltav) {
+    while (next_fdeltav > fdeltav && next_fdeltav > epsilon) {
       float g0 = sqrtf(fdeltav), g1 = sqrtf(next_fdeltav);
       float g01 = -g0 / len_v3(step);
       float det = 2.0f * (g1 - g0 - g01);