Math Lib: optimize axis_dominant_v3_to_m3, approx 6x speedup

build the matrix directly rather then calculating with axis/angle

also remove unused function calc_poly_plane

6 files changed, 57 insertions, 143 deletions

diff --git a/extern/carve/carve-util.cc b/extern/carve/carve-util.cc
index b268dae9dd6..d02b786fd2a 100644
--- a/extern/carve/carve-util.cc
+++ b/extern/carve/carve-util.cc
@@ -49,84 +49,53 @@ namespace {
 
 // Functions adopted from BLI_math.h to use Carve Vector and Matrix.
 
-void axis_angle_normalized_to_mat3(const Vector &normal,
-                                   const double angle,
-                                   Matrix3 *matrix)
+void transpose_m3__bli(double mat[3][3])
 {
-	double nsi[3], co, si, ico;
-
-	/* now convert this to a 3x3 matrix */
-	co = cos(angle);
-	si = sin(angle);
-
-	ico = (1.0 - co);
-	nsi[0] = normal[0] * si;
-	nsi[1] = normal[1] * si;
-	nsi[2] = normal[2] * si;
-
-	matrix->m[0][0] = ((normal[0] * normal[0]) * ico) + co;
-	matrix->m[0][1] = ((normal[0] * normal[1]) * ico) + nsi[2];
-	matrix->m[0][2] = ((normal[0] * normal[2]) * ico) - nsi[1];
-	matrix->m[1][0] = ((normal[0] * normal[1]) * ico) - nsi[2];
-	matrix->m[1][1] = ((normal[1] * normal[1]) * ico) + co;
-	matrix->m[1][2] = ((normal[1] * normal[2]) * ico) + nsi[0];
-	matrix->m[2][0] = ((normal[0] * normal[2]) * ico) + nsi[1];
-	matrix->m[2][1] = ((normal[1] * normal[2]) * ico) - nsi[0];
-	matrix->m[2][2] = ((normal[2] * normal[2]) * ico) + co;
+	double t;
+
+	t = mat[0][1];
+	mat[0][1] = mat[1][0];
+	mat[1][0] = t;
+	t = mat[0][2];
+	mat[0][2] = mat[2][0];
+	mat[2][0] = t;
+	t = mat[1][2];
+	mat[1][2] = mat[2][1];
+	mat[2][1] = t;
 }
 
-void axis_angle_to_mat3(const Vector &axis,
-                        const double angle,
-                        Matrix3 *matrix)
+void ortho_basis_v3v3_v3__bli(double r_n1[3], double r_n2[3], const double n[3])
 {
-	if (axis.length2() < FLT_EPSILON) {
-		*matrix = Matrix3();
-		return;
-	}
-
-	Vector nor = axis;
-	nor.normalize();
+	const double eps = FLT_EPSILON;
+	const double f = (n[0] * n[0]) + (n[1] * n[1]);
 
-	axis_angle_normalized_to_mat3(nor, angle, matrix);
-}
+	if (f > eps) {
+		const double d = 1.0f / sqrt(f);
 
-inline double saacos(double fac)
-{
-	if      (fac <= -1.0) return M_PI;
-	else if (fac >=  1.0) return 0.0;
-	else                  return acos(fac);
+		r_n1[0] =  n[1] * d;
+		r_n1[1] = -n[0] * d;
+		r_n1[2] =  0.0f;
+		r_n2[0] = -n[2] * r_n1[1];
+		r_n2[1] =  n[2] * r_n1[0];
+		r_n2[2] =  n[0] * r_n1[1] - n[1] * r_n1[0];
+	}
+	else {
+		/* degenerate case */
+		r_n1[0] = (n[2] < 0.0f) ? -1.0f : 1.0f;
+		r_n1[1] = r_n1[2] = r_n2[0] = r_n2[2] = 0.0f;
+		r_n2[1] = 1.0f;
+	}
 }
 
-bool axis_dominant_v3_to_m3(const Vector &normal,
-                            Matrix3 *matrix)
+void axis_dominant_v3_to_m3__bli(Matrix3 *r_mat, const Vector &normal)
 {
-	Vector up;
-	Vector axis;
-	double angle;
-
-	up.x = 0.0;
-	up.y = 0.0;
-	up.z = 1.0;
-
-	axis = carve::geom::cross(normal, up);
-	angle = saacos(carve::geom::dot(normal, up));
-
-	if (angle >= FLT_EPSILON) {
-		if (axis.length2() < FLT_EPSILON) {
-			axis[0] = 0.0;
-			axis[1] = 1.0;
-			axis[2] = 0.0;
-		}
+	memcpy(r_mat->m[2], normal.v, sizeof(double[3]));
+	ortho_basis_v3v3_v3__bli(r_mat->m[0], r_mat->m[1], r_mat->m[2]);
 
-		axis_angle_to_mat3(axis, angle, matrix);
-		return true;
-	}
-	else {
-		*matrix = Matrix3();
-		return false;
-	}
+	transpose_m3__bli(r_mat->m);
 }
 
+
 void meshset_minmax(const MeshSet<3> *mesh,
                     Vector *min,
                     Vector *max)
@@ -581,7 +550,7 @@ bool carve_checkPolyPlanarAndGetNormal(const std::vector<Vector> &vertices,
 			double magnitude = normal.length2();
 
 			normal.normalize();
-			axis_dominant_v3_to_m3(normal, axis_matrix_r);
+			axis_dominant_v3_to_m3__bli(axis_matrix_r, normal);
 
 			Vector first_projected = *axis_matrix_r * vertices[verts_of_poly[0]];
 			double min_z = first_projected[2], max_z = first_projected[2];
diff --git a/source/blender/blenlib/BLI_math_geom.h b/source/blender/blenlib/BLI_math_geom.h
index 3c1b89fbd44..6ba8910db14 100644
--- a/source/blender/blenlib/BLI_math_geom.h
+++ b/source/blender/blenlib/BLI_math_geom.h
@@ -303,7 +303,7 @@ bool form_factor_visible_quad(const float p[3], const float n[3],
 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_to_m3(float r_mat[3][3], const float normal[3]);
 
 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;
diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c
index 45e8c425f64..c2e4799d7cb 100644
--- a/source/blender/blenlib/intern/math_geom.c
+++ b/source/blender/blenlib/intern/math_geom.c
@@ -2133,32 +2133,20 @@ void fill_poly_v2i_n(
  * \param r_mat The matrix to return.
  * \param normal A unit length vector.
  */
-bool axis_dominant_v3_to_m3(float r_mat[3][3], const float normal[3])
+void axis_dominant_v3_to_m3(float r_mat[3][3], const float normal[3])
 {
-	float up[3] = {0.0f, 0.0f, 1.0f};
-	float axis[3];
-	float angle;
-
-	/* double check they are normalized */
 	BLI_ASSERT_UNIT_V3(normal);
 
-	cross_v3_v3v3(axis, normal, up);
-	angle = saacos(dot_v3v3(normal, up));
+	copy_v3_v3(r_mat[2], normal);
+	ortho_basis_v3v3_v3(r_mat[0], r_mat[1], r_mat[2]);
 
-	if (angle >= FLT_EPSILON) {
-		if (len_squared_v3(axis) < FLT_EPSILON) {
-			axis[0] = 0.0f;
-			axis[1] = 1.0f;
-			axis[2] = 0.0f;
-		}
+	BLI_ASSERT_UNIT_V3(r_mat[0]);
+	BLI_ASSERT_UNIT_V3(r_mat[1]);
 
-		axis_angle_to_mat3(r_mat, axis, angle);
-		return true;
-	}
-	else {
-		unit_m3(r_mat);
-		return false;
-	}
+	transpose_m3(r_mat);
+
+	BLI_assert(!is_negative_m3(r_mat));
+	BLI_assert(fabsf(dot_m3_v3_row_z(r_mat, normal) - 1.0f) < BLI_ASSERT_UNIT_EPSILON);
 }
 
 /****************************** Interpolation ********************************/
diff --git a/source/blender/bmesh/intern/bmesh_interp.c b/source/blender/bmesh/intern/bmesh_interp.c
index e5921bb49d7..5bd1a72b5a2 100644
--- a/source/blender/bmesh/intern/bmesh_interp.c
+++ b/source/blender/bmesh/intern/bmesh_interp.c
@@ -317,12 +317,7 @@ static int quad_co(float *x, float *y, float v1[3], float v2[3], float v3[3], fl
 
 	/* rotate */
 	poly_rotate_plane(n, projverts, 5);
-	
-	/* flatten */
-	for (i = 0; i < 5; i++) {
-		projverts[i][2] = 0.0f;
-	}
-	
+
 	/* subtract origin */
 	for (i = 0; i < 4; i++) {
 		sub_v3_v3(projverts[i], projverts[4]);
diff --git a/source/blender/bmesh/intern/bmesh_polygon.c b/source/blender/bmesh/intern/bmesh_polygon.c
index 17996176555..a2cfdc41281 100644
--- a/source/blender/bmesh/intern/bmesh_polygon.c
+++ b/source/blender/bmesh/intern/bmesh_polygon.c
@@ -362,46 +362,6 @@ void BM_face_calc_center_mean_weighted(BMFace *f, float r_cent[3])
 }
 
 /**
- * COMPUTE POLY PLANE
- *
- * Projects a set polygon's vertices to
- * a plane defined by the average
- * of its edges cross products
- */
-void calc_poly_plane(float (*verts)[3], const int nverts)
-{
-	
-	float avgc[3], norm[3], mag, avgn[3];
-	float *v1, *v2, *v3;
-	int i;
-	
-	if (nverts < 3)
-		return;
-
-	zero_v3(avgn);
-	zero_v3(avgc);
-
-	for (i = 0; i < nverts; i++) {
-		v1 = verts[i];
-		v2 = verts[(i + 1) % nverts];
-		v3 = verts[(i + 2) % nverts];
-		normal_tri_v3(norm, v1, v2, v3);
-
-		add_v3_v3(avgn, norm);
-	}
-
-	if (UNLIKELY(normalize_v3(avgn) == 0.0f)) {
-		avgn[2] = 1.0f;
-	}
-	
-	for (i = 0; i < nverts; i++) {
-		v1 = verts[i];
-		mag = dot_v3v3(v1, avgn);
-		madd_v3_v3fl(v1, avgn, -mag);
-	}
-}
-
-/**
  * \brief BM LEGAL EDGES
  *
  * takes in a face and a list of edges, and sets to NULL any edge in
@@ -430,15 +390,18 @@ static void scale_edge_v2f(float v1[2], float v2[2], const float fac)
  * Rotates a polygon so that it's
  * normal is pointing towards the mesh Z axis
  */
-void poly_rotate_plane(const float normal[3], float (*verts)[3], const int nverts)
+void poly_rotate_plane(const float normal[3], float (*verts)[3], const unsigned int nverts)
 {
 	float mat[3][3];
+	float co[3];
+	unsigned int i;
 
-	if (axis_dominant_v3_to_m3(mat, normal)) {
-		int i;
-		for (i = 0; i < nverts; i++) {
-			mul_m3_v3(mat, verts[i]);
-		}
+	co[2] = 0.0f;
+
+	axis_dominant_v3_to_m3(mat, normal);
+	for (i = 0; i < nverts; i++) {
+		mul_v2_m3v3(co, mat, verts[i]);
+		copy_v3_v3(verts[i], co);
 	}
 }
 
diff --git a/source/blender/bmesh/intern/bmesh_private.h b/source/blender/bmesh/intern/bmesh_private.h
index 7ec418b2253..cac4713c8b2 100644
--- a/source/blender/bmesh/intern/bmesh_private.h
+++ b/source/blender/bmesh/intern/bmesh_private.h
@@ -72,8 +72,7 @@ enum {
 #define BM_ELEM_API_FLAG_TEST(element, f)    ((element)->head.api_flag &   (f))
 #define BM_ELEM_API_FLAG_CLEAR(element)      ((element)->head.api_flag = 0)
 
-void calc_poly_plane(float (*verts)[3], const int nverts);
-void poly_rotate_plane(const float normal[3], float (*verts)[3], const int nverts);
+void poly_rotate_plane(const float normal[3], float (*verts)[3], unsigned const int nverts);
 
 /* include the rest of our private declarations */
 #include "bmesh_structure.h"