Commit: 8191e65b58ff88e0edb6b8bc7908822f7465c206 Author: Alexander Gavrilov Date: Sat Sep 3 17:03:11 2022 +0300 Branches: temp-angavrilov https://developer.blender.org/rB8191e65b58ff88e0edb6b8bc7908822f7465c206
Shrinkwrap: fix stability of the Target Normal Project mode. This mode works by using an iterative process to solve a system of equations for each triangle to find a point on its surface that has the smooth normal pointing at the original point. If a point within the triangle is not found, the next triangle is searched. All instability with vertices jumping to the opposite side of the mesh is caused by incorrectly discarding triangles for various reasons when the solution is close to the triangle edge. In order to optimize performance the old code was aggressively aborting iteration when the local gradient at the edge was pointing outside domain. However, it is wrong because it can be caused by a sharp valley diagonal to the domain boundary with the bottom gently sloping towards a minimum within the domain. Now iteration is only aborted if the solution deviates nonsensically far from the domain. Otherwise, the iteration proceeds as usual, and the final result is checked against domain borders with epsilon. In addition, iterations can now be aborted based on the value of the Jacobian determinant, or the number of linear search correction steps. Both can signify a singularity, which can be caused by the lack of a real solution to the equation. Finally, custom correction clearly has to be done after the linear search phase of the iterative solver, because the linear correction math assumes running after a normal Newton method step, not some kind of custom clamping. Differential Revision: https://developer.blender.org/D15892 =================================================================== M source/blender/blenkernel/intern/shrinkwrap.cc M source/blender/blenlib/BLI_math_solvers.h M source/blender/blenlib/intern/math_solvers.c =================================================================== diff --git a/source/blender/blenkernel/intern/shrinkwrap.cc b/source/blender/blenkernel/intern/shrinkwrap.cc index e3ae489b966..db5614abaaf 100644 --- a/source/blender/blenkernel/intern/shrinkwrap.cc +++ b/source/blender/blenkernel/intern/shrinkwrap.cc @@ -776,6 +776,16 @@ static void target_project_tri_jacobian(void *userdata, const float x[3], float madd_v3_v3fl(r_jacobian[2], data->n1_minus_n2, x[1]); } +/* Retrieve maximum deviation of the barycentric weights outside the triangle. */ +static float target_project_tri_error(float x[3]) +{ + float error = 0.0f; + error = max_ff(error, -x[0]); + error = max_ff(error, -x[1]); + error = max_ff(error, x[0] + x[1] - 1.0f); + return error; +} + /* Clamp barycentric weights to the triangle. */ static void target_project_tri_clamp(float x[3]) { @@ -797,71 +807,26 @@ static bool target_project_tri_correct(void * /*userdata*/, float step[3], float x_next[3]) { - /* Insignificant correction threshold */ - 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; - - /* 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]); - - if (sum + sstep > 1.0f) { - float ldist = 1.0f - sum; - - /* If already at the boundary, slide along it. */ - if (ldist < epsilon * float(M_SQRT2)) { - float step_len = len_v2(step); - - /* Abort if the solution is clearly outside the domain. */ - 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, ldist / sstep); - fixed = true; - } - } - - /* 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; - } + const float error = target_project_tri_error(x_next); - /* Reset precision errors to stay at the boundary. */ - step[i] = x[i]; - fixed = true; - } - else { - /* Scale a significant step down to arrive at the boundary. */ - mul_v3_fl(step, x[i] / step[i]); - fixed = true; - } - } - } - - /* Recompute and clamp the new coordinates after step correction. */ - if (fixed) { - sub_v3_v3v3(x_next, x, step); - target_project_tri_clamp(x_next); + /* Immediately abort on a clearly wrong point. + * + * It is not appropriate to abort unless the value is extremely + * wrong, because if the goal function gradient forms a sharp + * valley with a gently sloping bottom, the iterative process may + * be dragged against the edge of the domain by local gradients + * for a number of steps even when the ultimate minimum is within + * the domain. + * + * This threshold basically represents an estimate of something so + * far outside the domain that all assumptions about the solution + * behavior are off. + */ + if (error >= 10.0f) { + return false; } + /* Otherwise just continue iterations normally. */ return true; } @@ -878,6 +843,7 @@ static bool target_project_solve_point_tri(const float *vtri_co[3], float magnitude_estimate = dist + len_manhattan_v3(vtri_co[0]) + len_manhattan_v3(vtri_co[1]) + len_manhattan_v3(vtri_co[2]); float epsilon = magnitude_estimate * 1.0e-6f; + float det_epsilon = square_f(magnitude_estimate) * 1e-6f; /* Initial solution vector: barycentric weights plus distance along normal. */ interp_weights_tri_v3(x, UNPACK3(vtri_co), hit_co); @@ -912,11 +878,14 @@ static bool target_project_solve_point_tri(const float *vtri_co[3], &tri_data, epsilon, 20, + det_epsilon, + 5, /* >= 32x overshoot */ trace, x, x); - if (ok) { + /* If found a suitable solution within the triangle. */ + if (ok && target_project_tri_error(x) < epsilon) { copy_v3_v3(r_hit_co, tri_data.co_interp); copy_v3_v3(r_hit_no, tri_data.no_interp); @@ -1093,7 +1062,7 @@ void BKE_shrinkwrap_find_nearest_surface(ShrinkwrapTreeData *tree, if (type == MOD_SHRINKWRAP_TARGET_PROJECT) { #ifdef TRACE_TARGET_PROJECT - printf("\n====== TARGET PROJECT START ======\n"); + printf("\n====== TARGET PROJECT START: (%.3f,%.3f,%.3f) ======\n", UNPACK3(co)); #endif BLI_bvhtree_find_nearest_ex( diff --git a/source/blender/blenlib/BLI_math_solvers.h b/source/blender/blenlib/BLI_math_solvers.h index ec86ee7f4e4..09e68ea3e52 100644 --- a/source/blender/blenlib/BLI_math_solvers.h +++ b/source/blender/blenlib/BLI_math_solvers.h @@ -91,6 +91,8 @@ typedef bool (*Newton3D_CorrectionFunc)(void *userdata, * \param userdata: Data for the callbacks. * \param epsilon: Desired precision. * \param max_iterations: Limit on the iterations. + * \param det_epsilon: Minimum allowed value for the Jacobian determinant (singularity threshold). + * \param max_line_steps: Maximum number of line search steps (overshoot by >= 2^steps times). * \param trace: Enables logging to console. * \param x_init: Initial solution vector. * \param result: Final result. @@ -102,6 +104,8 @@ bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta, void *userdata, float epsilon, int max_iterations, + float det_epsilon, + int max_line_steps, bool trace, const float x_init[3], float result[3]); diff --git a/source/blender/blenlib/intern/math_solvers.c b/source/blender/blenlib/intern/math_solvers.c index 1196c0bed7f..b2b79717367 100644 --- a/source/blender/blenlib/intern/math_solvers.c +++ b/source/blender/blenlib/intern/math_solvers.c @@ -153,12 +153,14 @@ bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta, void *userdata, float epsilon, int max_iterations, + float det_epsilon, + int max_line_steps, bool trace, const float x_init[3], float result[3]) { float fdelta[3], fdeltav, next_fdeltav; - float jacobian[3][3], step[3], x[3], x_next[3]; + float jacobian[3][3], inv_jacobian[3][3], step[3], x[3], x_next[3], x_check[3]; epsilon *= epsilon; @@ -175,24 +177,17 @@ bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta, /* Newton's method step. */ func_jacobian(userdata, x, jacobian); - if (!invert_m3(jacobian)) { + if (!invert_m3_m3_ex(inv_jacobian, jacobian, det_epsilon)) { + if (trace) { + /* The Jacobian determinant was close to zero, signifying a nearby singularity. */ + printf("SINGULARITY\n"); + } return false; } - mul_v3_m3v3(step, jacobian, fdelta); + mul_v3_m3v3(step, inv_jacobian, fdelta); sub_v3_v3v3(x_next, x, step); - /* Custom out-of-bounds value correction. */ - if (func_correction) { - if (trace) { - printf("%3d * (%g, %g, %g)\n", i, x_next[0], x_next[1], x_next[2]); - } - - if (!func_correction(userdata, x, step, x_next)) { - return false; - } - } - func_delta(userdata, x_next, fdelta); next_fdeltav = len_squared_v3(fdelta); @@ -201,7 +196,14 @@ bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta, } /* Line search correction. */ - while (next_fdeltav > fdeltav && next_fdeltav > epsilon) { + for (int j = 0; next_fdeltav > fdeltav && next_fdeltav > epsilon; j++) { + if (j >= max_line_steps) { + /* This means the iteration step was more than 2^j times the correct distance, + * which is a sign of a nearby singularity making the solution unstable. */ + printf("OVERSHOOT\n"); + return false; + } + float g0 = sqrtf(fdeltav), g1 = sqrtf(next_fdeltav); float g01 = -g0 / len_v3(step); float det = 2.0f * (g1 - g0 - g01); @@ -219,6 +221,24 @@ bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta, } } + /* Custom out-of-bounds value correction. */ + if (func_correction) { + copy_v3_v3(x_check, x_next); + + if (!func_correction(userdata, x, step, x_next)) { + return false; + } + + if (!equals_v3v3(x_check, x_next)) { + func_delta(userdata, x_next, fdelta); + next_fdeltav = len_squared_v3(fdelta); + + if (trace) { + printf("%3d * (%g, %g, %g) %g\n", i, x_next[0], x_next[1], x_next[2], next_fdeltav); + } + } + } + copy_v3_v3(x, x_next); fdeltav = next_fdeltav; } _______________________________________________ Bf-blender-cvs mailing list Bf-blender-cvs@blender.org List details, subscription details or unsubscribe: https://lists.blender.org/mailman/listinfo/bf-blender-cvs
