Fix T39210: Grid Fill is generating mesh that's inconsistent with selected edge loops

Issue was in BLI's rotation_between_vecs_to_quat(), which did not handled correctly cases where both vectors are colinear.

Patch by Campbell Barton and me.

Issue originaly tracked down by Yan Shi, many thanks!

6 files changed, 89 insertions, 30 deletions

diff --git a/source/blender/blenlib/BLI_math_geom.h b/source/blender/blenlib/BLI_math_geom.h
index 09ff3585312..8b5190b9672 100644
--- a/source/blender/blenlib/BLI_math_geom.h
+++ b/source/blender/blenlib/BLI_math_geom.h
@@ -304,8 +304,10 @@ float form_factor_hemi_poly(float p[3], float n[3],
                             float v1[3], float v2[3], float v3[3], float v4[3]);
 
 bool  axis_dominant_v3_to_m3(float r_mat[3][3], const float normal[3]);
-void  axis_dominant_v3(int *r_axis_a, int *r_axis_b, const float axis[3]);
-float axis_dominant_v3_max(int *r_axis_a, int *r_axis_b, const float axis[3]) ATTR_WARN_UNUSED_RESULT;
+
+MINLINE void  axis_dominant_v3(int *r_axis_a, int *r_axis_b, const float axis[3]);
+MINLINE float axis_dominant_v3_max(int *r_axis_a, int *r_axis_b, const float axis[3]) ATTR_WARN_UNUSED_RESULT;
+MINLINE int axis_dominant_v3_single(const float vec[3]);
 
 MINLINE int max_axis_v3(const float vec[3]);
 MINLINE int min_axis_v3(const float vec[3]);
diff --git a/source/blender/blenlib/BLI_math_vector.h b/source/blender/blenlib/BLI_math_vector.h
index 22ff9a3ce6e..aef06f77d72 100644
--- a/source/blender/blenlib/BLI_math_vector.h
+++ b/source/blender/blenlib/BLI_math_vector.h
@@ -253,6 +253,7 @@ void project_v3_v3v3(float r[3], const float p[3], const float n[3]);
 void project_v3_plane(float v[3], const float n[3], const float p[3]);
 void reflect_v3_v3v3(float r[3], const float v[3], const float n[3]);
 void ortho_basis_v3v3_v3(float r1[3], float r2[3], const float a[3]);
+void ortho_v3_v3(float p[3], const float v[3]);
 void bisect_v3_v3v3v3(float r[3], const float a[3], const float b[3], const float c[3]);
 void rotate_v3_v3v3fl(float v[3], const float p[3], const float axis[3], const float angle);
 void rotate_normalized_v3_v3v3fl(float v[3], const float p[3], const float axis[3], const float angle);
diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c
index 2a19584cdb3..b1fd0fe5663 100644
--- a/source/blender/blenlib/intern/math_geom.c
+++ b/source/blender/blenlib/intern/math_geom.c
@@ -2131,31 +2131,6 @@ bool axis_dominant_v3_to_m3(float r_mat[3][3], const float normal[3])
 	}
 }
 
-/* get the 2 dominant axis values, 0==X, 1==Y, 2==Z */
-void axis_dominant_v3(int *r_axis_a, int *r_axis_b, const float axis[3])
-{
-	const float xn = fabsf(axis[0]);
-	const float yn = fabsf(axis[1]);
-	const float zn = fabsf(axis[2]);
-
-	if      (zn >= xn && zn >= yn) { *r_axis_a = 0; *r_axis_b = 1; }
-	else if (yn >= xn && yn >= zn) { *r_axis_a = 0; *r_axis_b = 2; }
-	else                           { *r_axis_a = 1; *r_axis_b = 2; }
-}
-
-/* same as axis_dominant_v3 but return the max value */
-float axis_dominant_v3_max(int *r_axis_a, int *r_axis_b, const float axis[3])
-{
-	const float xn = fabsf(axis[0]);
-	const float yn = fabsf(axis[1]);
-	const float zn = fabsf(axis[2]);
-
-	if      (zn >= xn && zn >= yn) { *r_axis_a = 0; *r_axis_b = 1; return zn; }
-	else if (yn >= xn && yn >= zn) { *r_axis_a = 0; *r_axis_b = 2; return yn; }
-	else                           { *r_axis_a = 1; *r_axis_b = 2; return xn; }
-}
-
-
 /****************************** Interpolation ********************************/
 
 static float tri_signed_area(const float v1[3], const float v2[3], const float v3[3], const int i, const int j)
diff --git a/source/blender/blenlib/intern/math_geom_inline.c b/source/blender/blenlib/intern/math_geom_inline.c
index a5906dbdf6c..c6f36464d51 100644
--- a/source/blender/blenlib/intern/math_geom_inline.c
+++ b/source/blender/blenlib/intern/math_geom_inline.c
@@ -155,6 +155,41 @@ MINLINE void madd_sh_shfl(float r[9], const float sh[9], const float f)
 	add_sh_shsh(r, r, tmp);
 }
 
+/* get the 2 dominant axis values, 0==X, 1==Y, 2==Z */
+MINLINE void axis_dominant_v3(int *r_axis_a, int *r_axis_b, const float axis[3])
+{
+	const float xn = fabsf(axis[0]);
+	const float yn = fabsf(axis[1]);
+	const float zn = fabsf(axis[2]);
+
+	if      (zn >= xn && zn >= yn) { *r_axis_a = 0; *r_axis_b = 1; }
+	else if (yn >= xn && yn >= zn) { *r_axis_a = 0; *r_axis_b = 2; }
+	else                           { *r_axis_a = 1; *r_axis_b = 2; }
+}
+
+/* same as axis_dominant_v3 but return the max value */
+MINLINE float axis_dominant_v3_max(int *r_axis_a, int *r_axis_b, const float axis[3])
+{
+	const float xn = fabsf(axis[0]);
+	const float yn = fabsf(axis[1]);
+	const float zn = fabsf(axis[2]);
+
+	if      (zn >= xn && zn >= yn) { *r_axis_a = 0; *r_axis_b = 1; return zn; }
+	else if (yn >= xn && yn >= zn) { *r_axis_a = 0; *r_axis_b = 2; return yn; }
+	else                           { *r_axis_a = 1; *r_axis_b = 2; return xn; }
+}
+
+/* get the single dominant axis value, 0==X, 1==Y, 2==Z */
+MINLINE int axis_dominant_v3_single(const float vec[3])
+{
+	const float x = fabsf(vec[0]);
+	const float y = fabsf(vec[1]);
+	const float z = fabsf(vec[2]);
+	return ((x > y) ?
+	       ((x > z) ? 0 : 2) :
+	       ((y > z) ? 1 : 2));
+}
+
 MINLINE int max_axis_v3(const float vec[3])
 {
 	const float x = vec[0];
diff --git a/source/blender/blenlib/intern/math_rotation.c b/source/blender/blenlib/intern/math_rotation.c
index c7c0626019f..d1753a77f7f 100644
--- a/source/blender/blenlib/intern/math_rotation.c
+++ b/source/blender/blenlib/intern/math_rotation.c
@@ -404,13 +404,29 @@ float normalize_qt_qt(float r[4], const float q[4])
 void rotation_between_vecs_to_quat(float q[4], const float v1[3], const float v2[3])
 {
 	float axis[3];
-	float angle;
 
 	cross_v3_v3v3(axis, v1, v2);
 
-	angle = angle_normalized_v3v3(v1, v2);
+	if (normalize_v3(axis) > FLT_EPSILON) {
+		float angle;
 
-	axis_angle_to_quat(q, axis, angle);
+		angle = angle_normalized_v3v3(v1, v2);
+
+		axis_angle_normalized_to_quat(q, axis, angle);
+	}
+	else {
+		/* degenerate case */
+
+		if (dot_v3v3(v1, v2) > 0.0f) {
+			/* Same vectors, zero rotation... */
+			unit_qt(q);
+		}
+		else {
+			/* Colinear but opposed vectors, 180 rotation... */
+			ortho_v3_v3(axis, v1);
+			axis_angle_to_quat(q, axis, (float)M_PI);
+		}
+	}
 }
 
 void rotation_between_quats_to_quat(float q[4], const float q1[4], const float q2[4])
diff --git a/source/blender/blenlib/intern/math_vector.c b/source/blender/blenlib/intern/math_vector.c
index 098272c9bc0..b189fc32284 100644
--- a/source/blender/blenlib/intern/math_vector.c
+++ b/source/blender/blenlib/intern/math_vector.c
@@ -514,6 +514,36 @@ void ortho_basis_v3v3_v3(float v1[3], float v2[3], const float v[3])
 	}
 }
 
+/**
+ * Calculates \a p - a perpendicular vector to \a v
+ *
+ * \note return vector won't maintain same length.
+ */
+void ortho_v3_v3(float p[3], const float v[3])
+{
+	const int axis = axis_dominant_v3_single(v);
+
+	BLI_assert(p != v);
+
+	switch (axis) {
+		case 0:
+			p[0] = -v[1] - v[2];
+			p[1] =  v[0];
+			p[2] =  v[0];
+			break;
+		case 1:
+			p[0] =  v[1];
+			p[1] = -v[0] - v[2];
+			p[2] =  v[1];
+			break;
+		case 2:
+			p[0] =  v[2];
+			p[1] =  v[2];
+			p[2] = -v[0] - v[1];
+			break;
+	}
+}
+
 /* Rotate a point p by angle theta around an arbitrary axis r
  * http://local.wasp.uwa.edu.au/~pbourke/geometry/
  */