diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h new file mode 100644 index 0000000000000000000000000000000000000000..4028729b4169383b4c45e561b195815d0a75e8f4 --- /dev/null +++ b/GeoLib/AnalyticalGeometry-impl.h @@ -0,0 +1,106 @@ +/** + * \copyright + * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +namespace GeoLib +{ + +template <class T_POINT> +void getNewellPlane (const std::vector<T_POINT*>& pnts, + MathLib::Vector3 &plane_normal, + double& d) +{ + d = 0; + MathLib::Vector3 centroid; + std::size_t n_pnts(pnts.size()); + for (size_t i(n_pnts - 1), j(0); j < n_pnts; i = j, j++) { + plane_normal[0] += ((*(pnts[i]))[1] - (*(pnts[j]))[1]) + * ((*(pnts[i]))[2] + (*(pnts[j]))[2]); // projection on yz + plane_normal[1] += ((*(pnts[i]))[2] - (*(pnts[j]))[2]) + * ((*(pnts[i]))[0] + (*(pnts[j]))[0]); // projection on xz + plane_normal[2] += ((*(pnts[i]))[0] - (*(pnts[j]))[0]) + * ((*(pnts[i]))[1] + (*(pnts[j]))[1]); // projection on xy + + centroid += *(pnts[j]); + } + + plane_normal.normalize(); + d = MathLib::scalarProduct(centroid, plane_normal) / n_pnts; +} + +template<class T_MATRIX> +void compute2DRotationMatrixToX(MathLib::Vector3 const& v, + T_MATRIX & rot_mat) +{ + const double cos_theta = v[0]; + const double sin_theta = v[1]; + rot_mat(0,0) = rot_mat(1,1) = cos_theta; + rot_mat(0,1) = sin_theta; + rot_mat(1,0) = -sin_theta; +} + +template <class T_MATRIX> +void compute3DRotationMatrixToX(MathLib::Vector3 const& v, + T_MATRIX & rot_mat) +{ + // a vector on the plane + MathLib::Vector3 yy(0, 0, 0); + if (fabs(v[0]) > 0.0 && fabs(v[1]) + fabs(v[2]) < std::numeric_limits<double>::epsilon()) { + yy[2] = 1.0; + } else if (fabs(v[1]) > 0.0 && fabs(v[0]) + fabs(v[2]) < std::numeric_limits<double>::epsilon()) { + yy[0] = 1.0; + } else if (fabs(v[2]) > 0.0 && fabs(v[0]) + fabs(v[1]) < std::numeric_limits<double>::epsilon()) { + yy[1] = 1.0; + } else { + for (unsigned i = 0; i < 3; i++) { + if (fabs(v[i]) > 0.0) { + yy[i] = -v[i]; + break; + } + } + } + // z"_vec + MathLib::Vector3 zz(MathLib::crossProduct(v, yy)); + zz.normalize(); + // y"_vec + yy = MathLib::crossProduct(zz, v); + yy.normalize(); + + for (unsigned i=0; i<3; ++i) { + rot_mat(0, i) = v[i]; + rot_mat(1, i) = yy[i]; + rot_mat(2, i) = zz[i]; + } +} + +template <class T_MATRIX> +void computeRotationMatrixToXY2(const std::vector<MathLib::Point3d*> &vec_pt, + T_MATRIX & rot_mat) +{ + assert(vec_pt.size()>2); + // x"_vec + MathLib::Vector3 xx((*vec_pt[0]), (*vec_pt[1])); + xx.normalize(); + // a vector on the plane + MathLib::Vector3 yy((*vec_pt[1]), (*vec_pt[2])); + // z"_vec + MathLib::Vector3 zz(MathLib::crossProduct(xx, yy)); + zz.normalize(); + // y"_vec + yy = MathLib::crossProduct(zz, xx); + yy.normalize(); + + for (size_t i=0; i<3; ++i) { + rot_mat(0, i) = xx[i]; + rot_mat(1, i) = yy[i]; + rot_mat(2, i) = zz[i]; + } +} + +} // end namespace GeoLib + diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp index d651df8b41aaab8f78dcbfabee6853333f3ac36a..bd9900cdd2e7c9983d0219a653172b873d471d71 100644 --- a/GeoLib/AnalyticalGeometry.cpp +++ b/GeoLib/AnalyticalGeometry.cpp @@ -348,27 +348,6 @@ double calcTetrahedronVolume(MathLib::Point3d const& x1, return std::abs(GeoLib::scalarTriple(ac, ad, ab)) / 6.0; } -// NewellPlane from book Real-Time Collision detection p. 494 -void getNewellPlane(const std::vector<GeoLib::Point*>& pnts, MathLib::Vector3 &plane_normal, double& d) -{ - d = 0; - MathLib::Vector3 centroid; - size_t n_pnts(pnts.size()); - for (size_t i(n_pnts - 1), j(0); j < n_pnts; i = j, j++) { - plane_normal[0] += ((*(pnts[i]))[1] - (*(pnts[j]))[1]) - * ((*(pnts[i]))[2] + (*(pnts[j]))[2]); // projection on yz - plane_normal[1] += ((*(pnts[i]))[2] - (*(pnts[j]))[2]) - * ((*(pnts[i]))[0] + (*(pnts[j]))[0]); // projection on xz - plane_normal[2] += ((*(pnts[i]))[0] - (*(pnts[j]))[0]) - * ((*(pnts[i]))[1] + (*(pnts[j]))[1]); // projection on xy - - centroid += *(pnts[j]); - } - - plane_normal *= 1.0 / plane_normal.getLength(); - d = MathLib::scalarProduct(centroid, plane_normal) / n_pnts; -} - void computeRotationMatrixToXY(MathLib::Vector3 const& n, MathLib::DenseMatrix<double> & rot_mat) { diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h index 1198202138b4b9b7464f52ccb4a3f7814a7f75a8..d3be61c4469be5d60a11cc64a66d5dce71480d14 100644 --- a/GeoLib/AnalyticalGeometry.h +++ b/GeoLib/AnalyticalGeometry.h @@ -60,10 +60,27 @@ Orientation getOrientation (const GeoLib::Point* p0, * @param plane_normal the normal of the plane the polygon is located in * @param d parameter from the plane equation */ -void getNewellPlane (const std::vector<GeoLib::Point*>& pnts, +template <class T_POINT> +void getNewellPlane (const std::vector<T_POINT*>& pnts, MathLib::Vector3 &plane_normal, double& d); +/** + * Computes a rotation matrix that rotates the given 2D normal vector parallel to X-axis + * @param v a 2D normal vector to be rotated + * @param rot_mat a 2x2 rotation matrix + */ +template<class T_MATRIX> +void compute2DRotationMatrixToX(MathLib::Vector3 const& v, T_MATRIX & rot_mat); + +/** + * Computes a rotation matrix that rotates the given 3D normal vector parallel to X-axis + * @param v a 3D normal vector to be rotated + * @param rot_mat a 3x3 rotation matrix + */ +template <class T_MATRIX> +void compute3DRotationMatrixToX(MathLib::Vector3 const& v, T_MATRIX & rot_mat); + /** * 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 @@ -72,6 +89,15 @@ void getNewellPlane (const std::vector<GeoLib::Point*>& pnts, 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$z\f$ axis. + * @param vec_pt a vector of 3d points that is rotated parallel to the \f$z\f$ axis + * @param rot_mat 3x3 rotation matrix + */ +template <class T_MATRIX> +void computeRotationMatrixToXY2(const std::vector<MathLib::Point3d*> &vec_pt, + T_MATRIX & 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 @@ -296,4 +322,6 @@ GeoLib::Polygon rotatePolygonToXY(GeoLib::Polygon const& polygon_in, } // end namespace GeoLib +#include "AnalyticalGeometry-impl.h" + #endif /* ANALYTICAL_GEOMETRY_H_ */