diff --git a/MathLib/AnalyticalGeometry.cpp b/MathLib/AnalyticalGeometry.cpp
index 170ebeeef4f36a24f3bfcb9e4506f2435fee85b9..c50b0fc9d8e4fd161063c8c3ecb4c748dddfb860 100644
--- a/MathLib/AnalyticalGeometry.cpp
+++ b/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
diff --git a/MathLib/AnalyticalGeometry.h b/MathLib/AnalyticalGeometry.h
index 9a8874509b67ef8f83dec00f5da160b2ac63e7a2..20aa333834327e79c53e1b503a82d9314d898f5a 100644
--- a/MathLib/AnalyticalGeometry.h
+++ b/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);