[Bf-blender-cvs] [f2f3c96d5b0] temp-angavrilov: Shrinkwrap: fix stability of the Target Normal Project mode.
Alexander Gavrilov
noreply at git.blender.org
Wed Sep 7 10:11:27 CEST 2022
Commit: f2f3c96d5b0cc5c6e648fdc0a5f4f81a57ec10b7
Author: Alexander Gavrilov
Date: Sat Sep 3 17:03:11 2022 +0300
Branches: temp-angavrilov
https://developer.blender.org/rBf2f3c96d5b0cc5c6e648fdc0a5f4f81a57ec10b7
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 either if the solution deviates
nonsensically far from the domain, or if the gradient is proven
to slope towards the outside perpendicular to the boundary.
Until either condition is met, values are simply clamped to
the domain.
In addition, 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.c
M source/blender/blenlib/intern/math_solvers.c
===================================================================
diff --git a/source/blender/blenkernel/intern/shrinkwrap.c b/source/blender/blenkernel/intern/shrinkwrap.c
index 3b656cf60ba..88063769a61 100644
--- a/source/blender/blenkernel/intern/shrinkwrap.c
+++ b/source/blender/blenkernel/intern/shrinkwrap.c
@@ -779,18 +779,30 @@ static void target_project_tri_jacobian(void *userdata, const float x[3], float
}
/* Clamp barycentric weights to the triangle. */
-static void target_project_tri_clamp(float x[3])
+static float target_project_tri_clamp(float x[3])
{
+ float error = 0.0f;
+
if (x[0] < 0.0f) {
+ error = max_ff(error, -x[0]);
+
x[0] = 0.0f;
}
+
if (x[1] < 0.0f) {
+ error = max_ff(error, -x[1]);
+
x[1] = 0.0f;
}
+
if (x[0] + x[1] > 1.0f) {
+ error = max_ff(error, x[0] + x[1] - 1.0f);
+
x[0] = x[0] / (x[0] + x[1]);
x[1] = 1.0f - x[0];
}
+
+ return error;
}
/* Correct the Newton's method step to keep the coordinates within the triangle. */
@@ -799,71 +811,36 @@ static bool target_project_tri_correct(void *UNUSED(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;
- }
+ const float error = target_project_tri_clamp(x_next);
- /* 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;
- }
+ /* 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 start to break down.
+ */
+ if (error > 1.0f) {
+ return false;
}
- /* 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;
- }
+ if (error > 0) {
+ sub_v3_v3v3(step, x, 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;
- }
+ /* Abort if after clamping the step is reverted to effectively nothing,
+ * meaning the gradient is perpendicular to the boundary and the solution
+ * is likely truly outside the domain. */
+ if (error > 1e-4f && len_manhattan_v3(step) < 0.01f * error) {
+ return false;
}
}
- /* Recompute and clamp the new coordinates after step correction. */
- if (fixed) {
- sub_v3_v3v3(x_next, x, step);
- target_project_tri_clamp(x_next);
- }
-
return true;
}
@@ -1096,7 +1073,7 @@ void BKE_shrinkwrap_find_nearest_surface(struct 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/intern/math_solvers.c b/source/blender/blenlib/intern/math_solvers.c
index b5650410a70..af577f4fd71 100644
--- a/source/blender/blenlib/intern/math_solvers.c
+++ b/source/blender/blenlib/intern/math_solvers.c
@@ -158,7 +158,7 @@ bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta,
float result[3])
{
float fdelta[3], fdeltav, next_fdeltav;
- float jacobian[3][3], step[3], x[3], x_next[3];
+ float jacobian[3][3], step[3], x[3], x_next[3], x_check[3];
epsilon *= epsilon;
@@ -182,17 +182,6 @@ bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta,
mul_v3_m3v3(step, 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);
@@ -219,6 +208,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