diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp index 5e5bcb449139bffa4b0bd245707f895f24e2b6a8..37771bc670bbc9f4d9a0b1dc7d929204b01cdd22 100644 --- a/GeoLib/AnalyticalGeometry.cpp +++ b/GeoLib/AnalyticalGeometry.cpp @@ -269,58 +269,6 @@ void getNewellPlane(const std::vector<GeoLib::Point*>& pnts, MathLib::Vector3 &p d = MathLib::scalarProduct(centroid, plane_normal) / n_pnts; } -void rotatePointsToXY(MathLib::Vector3 &plane_normal, std::vector<GeoLib::Point*> &pnts) -{ - double small_value(std::numeric_limits<double>::epsilon()); - if (fabs(plane_normal[0]) < small_value && fabs(plane_normal[1]) < small_value) - return; - - MathLib::DenseMatrix<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(MathLib::Vector3 &n, std::vector<GeoLib::Point*> &pnts) -{ - double small_value(std::numeric_limits<double>::epsilon()); - 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])); - - MathLib::DenseMatrix<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(MathLib::Vector3 const& plane_normal, MathLib::DenseMatrix<double> & rot_mat) { // *** some frequently used terms *** @@ -344,6 +292,28 @@ void computeRotationMatrixToXY(MathLib::Vector3 const& plane_normal, MathLib::De rot_mat(2, 2) = plane_normal[2] * h2; } +void computeRotationMatrixToXZ(MathLib::Vector3 const& plane_normal, MathLib::DenseMatrix<double> & rot_mat) +{ + // *** some frequently used terms *** + // n_1^2 + n_2^2 + const double h0(plane_normal[0] * plane_normal[0] + plane_normal[1] * plane_normal[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 + plane_normal[2] * plane_normal[2])); + + // calc rotation matrix + rot_mat(0, 0) = plane_normal[1] * h1; + rot_mat(0, 1) = -plane_normal[0] * h1; + rot_mat(0, 2) = 0.0; + rot_mat(1, 0) = plane_normal[0] * h2; + rot_mat(1, 1) = plane_normal[1] * h2; + rot_mat(1, 2) = plane_normal[2] * h2; + rot_mat(2, 0) = plane_normal[0] * plane_normal[2] * h1 * h2; + rot_mat(2, 1) = plane_normal[1] * plane_normal[2] * h1 * h2; + rot_mat(2, 2) = -sqrt(h0) * h2; +} + void rotatePoints(MathLib::DenseMatrix<double> const& rot_mat, std::vector<GeoLib::Point*> &pnts) { double* tmp (NULL); @@ -356,6 +326,48 @@ void rotatePoints(MathLib::DenseMatrix<double> const& rot_mat, std::vector<GeoLi } } +void rotatePointsToXY(std::vector<GeoLib::Point*> &pnts) +{ + assert(pnts.size()>2); + // calculate supporting plane + MathLib::Vector3 plane_normal; + double d; + // compute the plane normal + GeoLib::getNewellPlane(pnts, plane_normal, d); + + const double tol (std::numeric_limits<double>::epsilon()); + if (std::abs(plane_normal[0]) > tol || std::abs(plane_normal[1]) > tol) { + // rotate copied points into x-y-plane + MathLib::DenseMatrix<double> rot_mat(3, 3); + computeRotationMatrixToXY(plane_normal, rot_mat); + rotatePoints(rot_mat, pnts); + } + + for (std::size_t k(0); k<pnts.size(); k++) + (*(pnts[k]))[2] = 0.0; // should be -= d but there are numerical errors +} + +void rotatePointsToXZ(std::vector<GeoLib::Point*> &pnts) +{ + assert(pnts.size()>2); + // calculate supporting plane + MathLib::Vector3 plane_normal; + double d; + // compute the plane normal + GeoLib::getNewellPlane(pnts, plane_normal, d); + + const double tol (std::numeric_limits<double>::epsilon()); + if (std::abs(plane_normal[0]) > tol || std::abs(plane_normal[1]) > tol) { + // rotate copied points into x-z-plane + MathLib::DenseMatrix<double> rot_mat(3, 3); + computeRotationMatrixToXZ(plane_normal, rot_mat); + rotatePoints(rot_mat, pnts); + } + + for (std::size_t k(0); k<pnts.size(); k++) + (*(pnts[k]))[1] = 0.0; // should be -= d but there are numerical errors +} + GeoLib::Point* triangleLineIntersection(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c, GeoLib::Point const& p, GeoLib::Point const& q) { const MathLib::Vector3 pq(p, q); diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h index 1635a10cdf467b23d07b9a5314bb35a0b33fcbb9..9aed682d823bf508076880a638387cc83735f726 100644 --- a/GeoLib/AnalyticalGeometry.h +++ b/GeoLib/AnalyticalGeometry.h @@ -59,26 +59,6 @@ void getNewellPlane (const std::vector<GeoLib::Point*>& pnts, MathLib::Vector3 &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::Vector3 &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. - * @sa getNewellPlane() - */ -void rotatePointsToXZ(MathLib::Vector3 &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 plane_normal the (3d) vector that is rotated parallel to the \f$z\f$ axis @@ -87,6 +67,14 @@ void rotatePointsToXZ(MathLib::Vector3 &plane_normal, std::vector<GeoLib::Point* void computeRotationMatrixToXY(MathLib::Vector3 const& plane_normal, MathLib::DenseMatrix<double> & rot_mat); +/** + * Method computes the rotation matrix that rotates the given vector parallel to the \f$y\f$ axis. + * @param plane_normal the (3d) vector that is rotated parallel to the \f$y\f$ axis + * @param rot_mat 3x3 rotation matrix + */ +void computeRotationMatrixToXZ(MathLib::Vector3 const& plane_normal, + MathLib::DenseMatrix<double> & rot_mat); + /** * rotate points according to the rotation matrix * @param rot_mat 3x3 dimensional rotation matrix @@ -94,6 +82,22 @@ void computeRotationMatrixToXY(MathLib::Vector3 const& plane_normal, */ void rotatePoints(MathLib::DenseMatrix<double> const& rot_mat, std::vector<GeoLib::Point*> &pnts); +/** + * rotate points to X-Y plane + * @param pnts a vector of points with a minimum length of three. + * Points are rotated using a rotation matrix computed from the first three points + * in the vector. Point coordinates are modified as a result of the rotation. + */ +void rotatePointsToXY(std::vector<GeoLib::Point*> &pnts); + +/** + * rotate points to X-Z plane + * @param pnts a vector of points with a minimum length of three. + * Points are rotated using a rotation matrix computed from the first three points + * in the vector. Point coordinates are modified as a result of the rotation. + */ +void rotatePointsToXZ(std::vector<GeoLib::Point*> &pnts); + bool isPointInTriangle (const GeoLib::Point* p, const GeoLib::Point* a, const GeoLib::Point* b, const GeoLib::Point* c); diff --git a/GeoLib/EarClippingTriangulation.cpp b/GeoLib/EarClippingTriangulation.cpp index c575423c50f7a65aa79ba8c0bb290b4010966871..80a4ec541378f745d582ce3bb30b2a96452718f1 100644 --- a/GeoLib/EarClippingTriangulation.cpp +++ b/GeoLib/EarClippingTriangulation.cpp @@ -32,7 +32,7 @@ EarClippingTriangulation::EarClippingTriangulation(const GeoLib::Polygon* polygo copyPolygonPoints (polygon); if (rot) { - rotate (); + rotatePointsToXY (_pnts); ensureCWOrientation (); } @@ -79,24 +79,6 @@ void EarClippingTriangulation::copyPolygonPoints (const GeoLib::Polygon* polygon } } -void EarClippingTriangulation::rotate () -{ - // calculate supporting plane - MathLib::Vector3 plane_normal; - double d; - // compute the plane normal - getNewellPlane(_pnts, plane_normal, d); - - double tol (std::numeric_limits<double>::epsilon()); - if (fabs(plane_normal[0]) > tol || fabs(plane_normal[1]) > tol) { - // rotate copied points into x-y-plane - rotatePointsToXY(plane_normal, _pnts); - } - - for (std::size_t k(0); k<_pnts.size(); k++) - (*(_pnts[k]))[2] = 0.0; // should be -= d but there are numerical errors -} - void EarClippingTriangulation::ensureCWOrientation () { std::size_t n_pnts (_pnts.size()); diff --git a/GeoLib/EarClippingTriangulation.h b/GeoLib/EarClippingTriangulation.h index a8dbcdb663523bf3edf0cc9629d64d7ed9df1f3e..270a5b1ca861ff69722c4e92c69df58c594ecc10 100644 --- a/GeoLib/EarClippingTriangulation.h +++ b/GeoLib/EarClippingTriangulation.h @@ -39,7 +39,6 @@ private: * copies the points of the polygon to the vector _pnts */ inline void copyPolygonPoints (const GeoLib::Polygon* polygon); - inline void rotate (); inline void ensureCWOrientation (); inline bool isEar(std::size_t v0, std::size_t v1, std::size_t v2) const; diff --git a/GeoLib/Polygon.cpp b/GeoLib/Polygon.cpp index a0ef14c5275be917a2847937a068b2596030fad8..93c4c394ea84f9df49d1fe0558c4e69a9bd3c3e5 100644 --- a/GeoLib/Polygon.cpp +++ b/GeoLib/Polygon.cpp @@ -241,16 +241,8 @@ void Polygon::ensureCWOrientation () for (std::size_t k(0); k < n_pnts; k++) tmp_polygon_pnts.push_back (new GeoLib::Point (*(this->getPoint(k)))); - // *** calculate supporting plane (plane normal and - MathLib::Vector3 plane_normal; - double d; - GeoLib::getNewellPlane(tmp_polygon_pnts, plane_normal, d); - - // *** rotate if necessary - double tol (std::numeric_limits<double>::epsilon()); - if (fabs(plane_normal[0]) > tol || fabs(plane_normal[1]) > tol) - // rotate copied points into x-y-plane - GeoLib::rotatePointsToXY(plane_normal, tmp_polygon_pnts); + // rotate copied points into x-y-plane + GeoLib::rotatePointsToXY(tmp_polygon_pnts); for (std::size_t k(0); k < tmp_polygon_pnts.size(); k++) (*(tmp_polygon_pnts[k]))[2] = 0.0; // should be -= d but there are numerical errors