diff options
author | Alexander Gavrilov <angavrilov@gmail.com> | 2019-10-20 00:15:10 +0300 |
---|---|---|
committer | Alexander Gavrilov <angavrilov@gmail.com> | 2019-10-20 16:53:21 +0300 |
commit | 47d7f5e20041707b18f021801517db970931fec2 (patch) | |
tree | 1f61f9bc8face4ff0c08be4ca09131c7d707e300 | |
parent | fed27c25aa146efa530f3172cf8318c9b468e3e2 (diff) |
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.
-rw-r--r-- | source/blender/blenkernel/intern/shrinkwrap.c | 48 | ||||
-rw-r--r-- | source/blender/blenlib/intern/math_solvers.c | 6 |
2 files changed, 36 insertions, 18 deletions
diff --git a/source/blender/blenkernel/intern/shrinkwrap.c b/source/blender/blenkernel/intern/shrinkwrap.c index 28f552cec2e..6f755aa6460 100644 --- a/source/blender/blenkernel/intern/shrinkwrap.c +++ b/source/blender/blenkernel/intern/shrinkwrap.c @@ -866,9 +866,9 @@ static bool target_project_solve_point_tri(const float *vtri_co[3], { float x[3], tmp[3]; float dist = sqrtf(hit_dist_sq); - float epsilon = dist * 1.0e-5f; - - CLAMP_MIN(epsilon, 1.0e-5f); + 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; /* Initial solution vector: barycentric weights plus distance along normal. */ interp_weights_tri_v3(x, UNPACK3(vtri_co), hit_co); @@ -929,7 +929,7 @@ static bool update_hit(BVHTreeNearest *nearest, if (dist_sq < nearest->dist_sq) { #ifdef TRACE_TARGET_PROJECT printf( - "===> %d (%.3f,%.3f,%.3f) %g < %g\n", index, UNPACK3(hit_co), dist_sq, nearest->dist_sq); + "#=#=#> %d (%.3f,%.3f,%.3f) %g < %g\n", index, UNPACK3(hit_co), dist_sq, nearest->dist_sq); #endif nearest->index = index; nearest->dist_sq = dist_sq; @@ -1088,14 +1088,14 @@ void BKE_shrinkwrap_find_nearest_surface(struct ShrinkwrapTreeData *tree, if (type == MOD_SHRINKWRAP_TARGET_PROJECT) { #ifdef TRACE_TARGET_PROJECT - printf("====== TARGET PROJECT START ======\n"); + printf("\n====== TARGET PROJECT START ======\n"); #endif BLI_bvhtree_find_nearest_ex( tree->bvh, co, nearest, mesh_looptri_target_project, tree, BVH_NEAREST_OPTIMAL_ORDER); #ifdef TRACE_TARGET_PROJECT - printf("====== TARGET PROJECT END: %d %g ======\n", nearest->index, nearest->dist_sq); + printf("====== TARGET PROJECT END: %d %g ======\n\n", nearest->index, nearest->dist_sq); #endif if (nearest->index < 0) { @@ -1260,7 +1260,10 @@ static void shrinkwrap_snap_with_side(float r_point_co[3], float forcesign, bool forcesnap) { - float dist = len_v3v3(point_co, hit_co); + float delta[3]; + sub_v3_v3v3(delta, point_co, hit_co); + + float dist = len_v3(delta); /* If exactly on the surface, push out along normal */ if (dist < FLT_EPSILON) { @@ -1273,13 +1276,28 @@ static void shrinkwrap_snap_with_side(float r_point_co[3], } /* Move to the correct side if needed */ else { - float delta[3]; - sub_v3_v3v3(delta, point_co, hit_co); - float dsign = signf(dot_v3v3(delta, hit_no) * forcesign); + float dsign = signf(dot_v3v3(delta, hit_no)); + + if (forcesign == 0.0f) { + forcesign = dsign; + } /* If on the wrong side or too close, move to correct */ - if (forcesnap || dsign * dist < goal_dist) { - interp_v3_v3v3(r_point_co, point_co, hit_co, (dist - goal_dist * dsign) / dist); + if (forcesnap || dsign * dist * forcesign < goal_dist) { + mul_v3_fl(delta, dsign / dist); + + /* At very small distance, blend in the hit normal to stabilize math. */ + float dist_epsilon = (fabsf(goal_dist) + len_manhattan_v3(hit_co)) * 1e-4f; + + if (dist < dist_epsilon) { +#ifdef TRACE_TARGET_PROJECT + printf("zero_factor %g = %g / %g\n", dist / dist_epsilon, dist, dist_epsilon); +#endif + + interp_v3_v3v3(delta, hit_no, delta, dist / dist_epsilon); + } + + madd_v3_v3v3fl(r_point_co, hit_co, delta, goal_dist * forcesign); } else { copy_v3_v3(r_point_co, point_co); @@ -1304,13 +1322,13 @@ void BKE_shrinkwrap_snap_point_to_surface(const struct ShrinkwrapTreeData *tree, const float point_co[3], float r_point_co[3]) { - float dist, tmp[3]; + float tmp[3]; switch (mode) { /* Offsets along the line between point_co and hit_co. */ case MOD_SHRINKWRAP_ON_SURFACE: - if (goal_dist != 0 && (dist = len_v3v3(point_co, hit_co)) > FLT_EPSILON) { - interp_v3_v3v3(r_point_co, point_co, hit_co, (dist - goal_dist) / dist); + if (goal_dist != 0) { + shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, 0, true); } else { copy_v3_v3(r_point_co, hit_co); 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); |