[Bf-blender-cvs] [e92905beabd] temp-angavrilov: Shrinkwrap: fix stability of the Target Normal Project mode.
Alexander Gavrilov
noreply at git.blender.org
Sat Feb 4 21:41:01 CET 2023
Commit: e92905beabd66e4749e5596dc5f9df6380c5ec51
Author: Alexander Gavrilov
Date: Sat Sep 3 17:03:11 2022 +0300
Branches: temp-angavrilov
https://developer.blender.org/rBe92905beabd66e4749e5596dc5f9df6380c5ec51
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 1a0b33ca2b1..5f5f29a1b36 100644
--- a/source/blender/blenkernel/intern/shrinkwrap.cc
+++ b/source/blender/blenkernel/intern/shrinkwrap.cc
@@ -752,6 +752,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])
{
@@ -773,71 +783,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;
}
@@ -854,6 +819,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);
@@ -888,11 +854,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);
@@ -1069,7 +1038,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;
}
More information about the Bf-blender-cvs
mailing list