added function rotatePoints() that is applying a 3x3 rotation matrix to the...

added function rotatePoints() that is applying a 3x3 rotation matrix to the points given within the vector pnts
refactored rotatePointToXY()
added function rotatePointsToXZ()

MathLib/AnalyticalGeometry.cpp

@@ -172,49 +172,92 @@ void getNewellPlane(const std::vector<GeoLib::Point*>& pnts, Vector &plane_norma
 	d = centroid.Dot(plane_normal) / n_pnts;
 }
 
-
 void rotatePointsToXY(Vector &plane_normal,
-		std::vector<GeoLib::Point*> &pnts)
+                      std::vector<GeoLib::Point*> &pnts)
 {
 	double small_value (sqrt (std::numeric_limits<double>::min()));
 	if (fabs(plane_normal[0]) < small_value && fabs(plane_normal[1]) < small_value)
 		return;
 
+	Matrix<double> rot_mat(3, 3);
+	computeRotationMatrixToXY(plane_normal, rot_mat);
+	rotatePoints(rot_mat, pnts);
+
+	double* tmp (rot_mat * plane_normal.getCoords());
+	for (std::size_t j(0); j < 3; j++)
+		plane_normal[j] = tmp[j];
+
+	delete [] tmp;
+}
+
+void rotatePointsToXZ(Vector &n, std::vector<GeoLib::Point*> &pnts)
+{
+	double small_value (sqrt (std::numeric_limits<double>::min()));
+	if (fabs(n[0]) < small_value && fabs(n[1]) < small_value)
+		return;
+
+	// *** some frequently used terms ***
+	// n_1^2 + n_2^2
+	const double h0(n[0] * n[0] + n[1] * n[1]);
+	// 1 / sqrt (n_1^2 + n_2^2)
+	const double h1(1.0 / sqrt(h0));
+	// 1 / sqrt (n_1^2 + n_2^2 + n_3^2)
+	const double h2(1.0 / sqrt(h0 + n[2] * n[2]));
+
+	Matrix<double> rot_mat(3, 3);
+	// calc rotation matrix
+	rot_mat(0, 0) = n[1] * h1;
+	rot_mat(0, 1) = - n[0] * h1;
+	rot_mat(0, 2) = 0.0;
+	rot_mat(1, 0) = n[0] * h2;
+	rot_mat(1, 1) = n[1] * h2;
+	rot_mat(1, 2) = n[2] * h2;
+	rot_mat(2, 0) = n[0] * n[2] * h1 * h2;
+	rot_mat(2, 1) = n[1] * n[2] * h1 * h2;
+	rot_mat(2, 2) = - sqrt(h0) * h2;
+
+	rotatePoints(rot_mat, pnts);
+
+	double *tmp (rot_mat * n.getCoords());
+	for (std::size_t j(0); j < 3; j++)
+		n[j] = tmp[j];
+
+	delete [] tmp;
+}
+
+void computeRotationMatrixToXY(Vector const& plane_normal, Matrix<double> & rot_mat)
+{
 	// *** some frequently used terms ***
 	// sqrt (v_1^2 + v_2^2)
 	double h0(sqrt(plane_normal[0] * plane_normal[0] + plane_normal[1]
-			* plane_normal[1]));
+	               * plane_normal[1]));
 	// 1 / sqrt (v_1^2 + v_2^2)
 	double h1(1 / h0);
 	// 1 / sqrt (h0 + v_3^2)
 	double h2(1.0 / sqrt(h0 + plane_normal[2] * plane_normal[2]));
 
-	Matrix<double> rot_mat(3, 3);
-	// calc rotation matrix
+	// calculate entries of rotation matrix
 	rot_mat(0, 0) = plane_normal[2] * plane_normal[0] * h2 * h1;
 	rot_mat(0, 1) = plane_normal[2] * plane_normal[1] * h2 * h1;
-	rot_mat(0, 2) = - h0 * h2;
+	rot_mat(0, 2) = -h0 * h2;
 	rot_mat(1, 0) = -plane_normal[1] * h1;
-	rot_mat(1, 1) = plane_normal[0] * h1;;
+	rot_mat(1, 1) = plane_normal[0] * h1;
 	rot_mat(1, 2) = 0.0;
 	rot_mat(2, 0) = plane_normal[0] * h2;
 	rot_mat(2, 1) = plane_normal[1] * h2;
 	rot_mat(2, 2) = plane_normal[2] * h2;
+}
 
-	double *tmp (NULL);
-	size_t n_pnts(pnts.size());
-	for (size_t k(0); k < n_pnts; k++) {
+void rotatePoints(Matrix<double> const& rot_mat, std::vector<GeoLib::Point*> &pnts)
+{
+	double* tmp (NULL);
+	const std::size_t n_pnts(pnts.size());
+	for (std::size_t k(0); k < n_pnts; k++) {
 		tmp = rot_mat * pnts[k]->getCoords();
-		for (size_t j(0); j < 3; j++)
+		for (std::size_t j(0); j < 3; j++)
 			(*(pnts[k]))[j] = tmp[j];
 		delete [] tmp;
 	}
-
-	tmp = rot_mat * plane_normal.getCoords();
-	for (size_t j(0); j < 3; j++)
-		plane_normal[j] = tmp[j];
-
-	delete [] tmp;
 }
 
 } // end namespace MathLib

MathLib/AnalyticalGeometry.h

@@ -15,6 +15,7 @@
 
 // MathLib
 #include "Vector3.h"
+#include "LinAlg/Dense/Matrix.h"
 // GeoLib
 #include "Triangle.h"
 
@@ -55,11 +56,40 @@ Orientation getOrientation (const GeoLib::Point* p0, const GeoLib::Point* p1, co
 void getNewellPlane (const std::vector<GeoLib::Point*>& pnts, MathLib::Vector &plane_normal, double& d);
 
 /**
- *
+ * The vector plane_normal should be the surface normal of the plane surface described
+ * by the points within the vector pnts. See function getNewellPlane() to get the
+ * "plane normal" of a point set. The method rotates both the plane normal and
+ * the points. The plane normal is rotated such that it is parallel to the \f$z\f$ axis.
+ * @param plane_normal the normal of the plane
+ * @param pnts pointers to points in a vector that should be rotated
+ * @sa getNewellPlane()
+ */
+void rotatePointsToXY(MathLib::Vector &plane_normal, std::vector<GeoLib::Point*> &pnts);
+
+/**
+ * The vector plane_normal should be the surface normal of the plane surface described
+ * by the points within the vector pnts. See function getNewellPlane() to get the
+ * "plane normal" of a point set. The method rotates both the plane normal and
+ * the points. The plane normal is rotated such that it is parallel to the \f$y\f$ axis.
  * @param plane_normal
  * @param pnts
+ * @sa getNewellPlane()
  */
-void rotatePointsToXY(MathLib::Vector &plane_normal, std::vector<GeoLib::Point*> &pnts);
+void rotatePointsToXZ(MathLib::Vector &plane_normal, std::vector<GeoLib::Point*> &pnts);
+
+/**
+ * Method computes the rotation matrix that rotates the given vector parallel to the \f$z\f$ axis.
+ * @param vec the (3d) vector that is rotated parallel to the \f$z\f$ axis
+ * @param rot_mat 3x3 rotation matrix
+ */
+void computeRotationMatrixToXY(Vector const& plane_normal, Matrix<double> & rot_mat);
+
+/**
+ * rotate points according to the rotation matrix
+ * @param rot_mat 3x3 dimensional rotation matrix
+ * @param pnts vector of points
+ */
+void rotatePoints(Matrix<double> const& rot_mat, std::vector<GeoLib::Point*> &pnts);
 
 bool isPointInTriangle (const GeoLib::Point* p,
 		const GeoLib::Point* a, const GeoLib::Point* b, const GeoLib::Point* c);