Resolve the opposite vector ambiguity in Damped Track constraint.

Damped Track by specification attempts to arrive at the desired
direction via the shortest rotation. However with opposite vectors
there are infinitely many valid 180 degree rotations. Currently
it gives up and does nothing.

I think that it would be more reasonable to resolve the ambiguity
arbitrarily, so that Damped Track won't have a weird dead zone.
To make it more predictable I use a local axis.

In addition, the singularity area vicinity has some floating
point precision problems that result in significant jitter.
This applies workarounds for two causes of instability.

Differential Revision: https://developer.blender.org/D3530

3 files changed, 41 insertions, 2 deletions

diff --git a/source/blender/blenkernel/intern/constraint.c b/source/blender/blenkernel/intern/constraint.c
index 2a5a0cf9ae7..7d861658993 100644
--- a/source/blender/blenkernel/intern/constraint.c
+++ b/source/blender/blenkernel/intern/constraint.c
@@ -3709,7 +3709,7 @@ static void damptrack_evaluate(bConstraint *con, bConstraintOb *cob, ListBase *t
 		 *	- the min/max wrappers around (obvec . tarvec) result (stored temporarily in rangle)
 		 *	  are used to ensure that the smallest angle is chosen
 		 */
-		cross_v3_v3v3(raxis, obvec, tarvec);
+		cross_v3_v3v3_hi_prec(raxis, obvec, tarvec);
 
 		rangle = dot_v3v3(obvec, tarvec);
 		rangle = acosf(max_ff(-1.0f, min_ff(1.0f, rangle)));
@@ -3717,7 +3717,35 @@ static void damptrack_evaluate(bConstraint *con, bConstraintOb *cob, ListBase *t
 		/* construct rotation matrix from the axis-angle rotation found above
 		 *	- this call takes care to make sure that the axis provided is a unit vector first
 		 */
-		axis_angle_to_mat3(rmat, raxis, rangle);
+		float norm = normalize_v3(raxis);
+
+		if (norm < FLT_EPSILON) {
+			/* if dot product is nonzero, while cross is zero, we have two opposite vectors!
+			 *  - this is an ambiguity in the math that needs to be resolved arbitrarily,
+			 *    or there will be a case where damped track strangely does nothing
+			 *  - to do that, rotate around a different local axis
+			 */
+			float tmpvec[3];
+
+			if (fabsf(rangle) < M_PI - 0.01f) {
+				return;
+			}
+
+			rangle = M_PI;
+			copy_v3_v3(tmpvec, track_dir_vecs[(data->trackflag + 1) % 6]);
+			mul_mat3_m4_v3(cob->matrix, tmpvec);
+			cross_v3_v3v3(raxis, obvec, tmpvec);
+
+			if (normalize_v3(raxis) == 0.0f) {
+				return;
+			}
+		}
+		else if (norm < 0.1f) {
+			/* near 0 and Pi arcsin has way better precision than arccos */
+			rangle = (rangle > M_PI_2) ? M_PI - asinf(norm) : asinf(norm);
+		}
+
+		axis_angle_normalized_to_mat3(rmat, raxis, rangle);
 
 		/* rotate the owner in the way defined by this rotation matrix, then reapply the location since
 		 * we may have destroyed that in the process of multiplying the matrix
diff --git a/source/blender/blenlib/BLI_math_vector.h b/source/blender/blenlib/BLI_math_vector.h
index 3f603311530..39625346756 100644
--- a/source/blender/blenlib/BLI_math_vector.h
+++ b/source/blender/blenlib/BLI_math_vector.h
@@ -169,6 +169,7 @@ MINLINE double dot_v3db_v3fl(const double a[3], const float b[3]) ATTR_WARN_UNUS
 
 MINLINE float cross_v2v2(const float a[2], const float b[2]) ATTR_WARN_UNUSED_RESULT;
 MINLINE void cross_v3_v3v3(float r[3], const float a[3], const float b[3]);
+MINLINE void cross_v3_v3v3_hi_prec(float r[3], const float a[3], const float b[3]);
 
 MINLINE void add_newell_cross_v3_v3v3(float n[3], const float v_prev[3], const float v_curr[3]);
 
diff --git a/source/blender/blenlib/intern/math_vector_inline.c b/source/blender/blenlib/intern/math_vector_inline.c
index 08687a1ab47..715e2e65c96 100644
--- a/source/blender/blenlib/intern/math_vector_inline.c
+++ b/source/blender/blenlib/intern/math_vector_inline.c
@@ -753,6 +753,16 @@ MINLINE void cross_v3_v3v3(float r[3], const float a[3], const float b[3])
 	r[2] = a[0] * b[1] - a[1] * b[0];
 }
 
+/* cross product suffers from severe precision loss when vectors are
+ * nearly parallel or opposite; doing the computation in double helps a lot */
+MINLINE void cross_v3_v3v3_hi_prec(float r[3], const float a[3], const float b[3])
+{
+	BLI_assert(r != a && r != b);
+	r[0] = (float)((double)a[1] * (double)b[2] - (double)a[2] * (double)b[1]);
+	r[1] = (float)((double)a[2] * (double)b[0] - (double)a[0] * (double)b[2]);
+	r[2] = (float)((double)a[0] * (double)b[1] - (double)a[1] * (double)b[0]);
+}
+
 /* Newell's Method */
 /* excuse this fairly specific function,
  * its used for polygon normals all over the place