diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 362b8451d8cd9ea3eee7ea96e578c8ad479c1148..c38b98e3fa2e521bbb02f9d08f0057db0ecbf269 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -250,6 +250,14 @@ double calcTriangleArea(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::
return 0.5 * w.getLength();
}
+double calcTetrahedronVolume(const double* x1, const double* x2, const double* x3, const double* x4)
+{
+ const MathLib::Vector3 ab(x1, x2);
+ const MathLib::Vector3 ac(x1, x3);
+ const MathLib::Vector3 ad(x1, x4);
+ 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)
{
@@ -459,4 +467,5 @@ void computeAndInsertAllIntersectionPoints(GeoLib::PointVec &pnt_vec,
}
}
+
} // end namespace GeoLib
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 9a705b5d1e2828dd27a2c94107a0a7b99696bc3a..5ca4905023f3c0a25791b7ad28513ffe8d9b7161 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -106,6 +106,12 @@ void rotatePointsToXZ(std::vector<GeoLib::Point*> &pnts);
*/
double calcTriangleArea(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c);
+/**
+ * Calculates the volume of a tetrahedron.
+ * The formula is V=1/6*|a(b x c)| with a=x1->x2, b=x1->x3 and c=x1->x4.
+ */
+double calcTetrahedronVolume(const double* x1, const double* x2, const double* x3, const double* x4);
+
bool isPointInTriangle (const GeoLib::Point* p,
const GeoLib::Point* a, const GeoLib::Point* b, const GeoLib::Point* c);
@@ -194,6 +200,7 @@ void computeAndInsertAllIntersectionPoints(
GeoLib::PointVec &pnt_vec,
std::vector<GeoLib::Polyline*> & plys);
+
} // end namespace GeoLib
#endif /* ANALYTICAL_GEOMETRY_H_ */
diff --git a/GeoLib/MinimalBoundingSphere.cpp b/GeoLib/MinimalBoundingSphere.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..19774ae808ae6e24019b752daedccfda94037e33
--- /dev/null
+++ b/GeoLib/MinimalBoundingSphere.cpp
@@ -0,0 +1,192 @@
+/**
+ * \file Calculation of a minimum bounding sphere for a vector of points
+ * \author Karsten Rink
+ * \date 2014-07-11
+ * \brief Implementation of the MinimalBoundingSphere class.
+ *
+ * \copyright
+ * Copyright (c) 2012-2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "MinimalBoundingSphere.h"
+
+#include <ctime>
+
+#include "MathTools.h"
+#include "AnalyticalGeometry.h"
+
+namespace GeoLib {
+
+MinimalBoundingSphere::MinimalBoundingSphere()
+: _radius(-1), _center(std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max())
+{
+}
+
+MinimalBoundingSphere::MinimalBoundingSphere(GeoLib::Point const& p, double radius)
+: _radius(radius), _center(p)
+{
+}
+
+MinimalBoundingSphere::MinimalBoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q)
+: _radius(std::numeric_limits<double>::epsilon()), _center(p)
+{
+ MathLib::Vector3 const a(p, q);
+
+ if (a.getLength() > 0)
+ {
+ MathLib::Vector3 const o(0.5*a);
+ _radius = o.getLength() + std::numeric_limits<double>::epsilon();
+ _center = MathLib::Vector3(p) + o;
+ }
+}
+
+MinimalBoundingSphere::MinimalBoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q, GeoLib::Point const& r)
+{
+ MathLib::Vector3 const a(p,r);
+ MathLib::Vector3 const b(p,q);
+
+ MathLib::Vector3 const cross_ab(crossProduct(a,b));
+
+ if (cross_ab.getLength() > 0)
+ {
+ double const denom = 2.0 * scalarProduct(cross_ab,cross_ab);
+ MathLib::Vector3 const o = (scalarProduct(b,b) * crossProduct(cross_ab, a)
+ + scalarProduct(a,a) * crossProduct(b, cross_ab))
+ * (1.0 / denom);
+ _radius = o.getLength() + std::numeric_limits<double>::epsilon();
+ _center = MathLib::Vector3(p) + o;
+ }
+ else
+ {
+ MinimalBoundingSphere two_pnts_sphere;
+ if (a.getLength() > b.getLength())
+ two_pnts_sphere = MinimalBoundingSphere(p,r);
+ else
+ two_pnts_sphere = MinimalBoundingSphere(p,q);
+ _radius = two_pnts_sphere.getRadius();
+ _center = two_pnts_sphere.getCenter();
+ }
+}
+
+MinimalBoundingSphere::MinimalBoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q, GeoLib::Point const& r, GeoLib::Point const& s)
+{
+ MathLib::Vector3 const a(p, q);
+ MathLib::Vector3 const b(p, r);
+ MathLib::Vector3 const c(p, s);
+
+ if (!GeoLib::isCoplanar(p, q, r, s))
+ {
+ double const denom = 2.0 * GeoLib::scalarTriple(a,b,c);
+ MathLib::Vector3 const o = (scalarProduct(c,c) * crossProduct(a,b)
+ + scalarProduct(b,b) * crossProduct(c,a)
+ + scalarProduct(a,a) * crossProduct(b,c))
+ * (1.0 / denom);
+
+ _radius = o.getLength() + std::numeric_limits<double>::epsilon();
+ _center = MathLib::Vector3(p) + o;
+ }
+ else
+ {
+ MinimalBoundingSphere const pqr(p, q , r);
+ MinimalBoundingSphere const pqs(p, q , s);
+ MinimalBoundingSphere const prs(p, r , s);
+ MinimalBoundingSphere const qrs(q, r , s);
+ _radius = pqr.getRadius();
+ _center = pqr.getCenter();
+ if (_radius < pqs.getRadius())
+ {
+ _radius = pqs.getRadius();
+ _center = pqs.getCenter();
+ }
+ if (_radius < prs.getRadius())
+ {
+ _radius = prs.getRadius();
+ _center = prs.getCenter();
+ }
+ if (_radius < qrs.getRadius())
+ {
+ _radius = qrs.getRadius();
+ _center = qrs.getCenter();
+ }
+ }
+}
+
+MinimalBoundingSphere::MinimalBoundingSphere(std::vector<GeoLib::Point*> const& points)
+: _radius(-1), _center(0,0,0)
+{
+ std::vector<GeoLib::Point*> sphere_points(points);
+ MinimalBoundingSphere const bounding_sphere = recurseCalculation(sphere_points, 0, sphere_points.size(), 0);
+ _center = bounding_sphere.getCenter();
+ _radius = bounding_sphere.getRadius();
+}
+
+MinimalBoundingSphere MinimalBoundingSphere::recurseCalculation(std::vector<GeoLib::Point*> sphere_points, std::size_t start_idx, std::size_t length, std::size_t n_boundary_points)
+{
+ MinimalBoundingSphere sphere;
+ switch(n_boundary_points)
+ {
+ case 0:
+ sphere = MinimalBoundingSphere();
+ break;
+ case 1:
+ sphere = MinimalBoundingSphere(*sphere_points[start_idx-1]);
+ break;
+ case 2:
+ sphere = MinimalBoundingSphere(*sphere_points[start_idx-1], *sphere_points[start_idx-2]);
+ break;
+ case 3:
+ sphere = MinimalBoundingSphere(*sphere_points[start_idx-1], *sphere_points[start_idx-2], *sphere_points[start_idx-3]);
+ break;
+ case 4:
+ sphere = MinimalBoundingSphere(*sphere_points[start_idx-1], *sphere_points[start_idx-2], *sphere_points[start_idx-3], *sphere_points[start_idx-4]);
+ return sphere;
+ }
+
+ for(std::size_t i=0; i<length; ++i)
+ {
+ // current point is located outside of sphere
+ if (sphere.pointDistanceSquared(*sphere_points[start_idx+i]) > 0)
+ {
+ if (i>start_idx)
+ {
+ GeoLib::Point* tmp = sphere_points[start_idx+i];
+ std::copy(sphere_points.begin() + start_idx, sphere_points.begin() + (start_idx + i), sphere_points.begin() + (start_idx + 1));
+ sphere_points[start_idx] = tmp;
+ }
+ sphere = recurseCalculation(sphere_points, start_idx+1, i, n_boundary_points+1);
+ }
+ }
+ return sphere;
+}
+
+double MinimalBoundingSphere::pointDistanceSquared(GeoLib::Point const& pnt) const
+{
+ return MathLib::sqrDist(_center.getCoords(), pnt.getCoords())-(_radius*_radius);
+}
+
+std::vector<GeoLib::Point*>* MinimalBoundingSphere::getRandomSpherePoints(std::size_t n_points) const
+{
+ std::vector<GeoLib::Point*> *pnts = new std::vector<GeoLib::Point*>;
+ pnts->reserve(n_points);
+ srand ( static_cast<unsigned>(time(NULL)) );
+
+ for (std::size_t k(0); k<n_points; ++k)
+ {
+ MathLib::Vector3 vec (0,0,0);
+ double sum (0);
+ for (unsigned i=0; i<3; ++i)
+ {
+ vec[i] = (double)rand()-(RAND_MAX/2.0);
+ sum+=(vec[i]*vec[i]);
+ }
+ double const fac (_radius/sqrt(sum));
+ pnts->push_back(new GeoLib::Point(_center[0]+vec[0]*fac, _center[1]+vec[1]*fac, _center[2]+vec[2]*fac));
+ }
+ return pnts;
+}
+
+}
diff --git a/GeoLib/MinimalBoundingSphere.h b/GeoLib/MinimalBoundingSphere.h
new file mode 100644
index 0000000000000000000000000000000000000000..cf3f48665ce381a4a19422a7c1d95e31620bbd5f
--- /dev/null
+++ b/GeoLib/MinimalBoundingSphere.h
@@ -0,0 +1,84 @@
+/**
+ * \file Calculation of a minimum bounding sphere for a vector of points
+ * \author Karsten Rink
+ * \date 2014-07-11
+ * \brief Definition of the MinimalBoundingSphere class.
+ *
+ * \copyright
+ * Copyright (c) 2012-2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ *
+ */
+
+#ifndef MINIMALBOUNDINGSPHERE_H_
+#define MINIMALBOUNDINGSPHERE_H_
+
+#include <vector>
+
+#include "Vector3.h"
+#include "Point.h"
+
+namespace GeoLib
+{
+
+/**
+ * Calculated center and radius of a minimal bounding sphere for a given number of geometric points.
+ */
+class MinimalBoundingSphere
+{
+public:
+ /// Copy constructor
+ MinimalBoundingSphere(MinimalBoundingSphere const& sphere) = default;
+ /// Point-Sphere
+ MinimalBoundingSphere(GeoLib::Point const& p, double radius = std::numeric_limits<double>::epsilon());
+ /// Bounding sphere using two points
+ MinimalBoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q);
+ /// Bounding sphere using three points
+ MinimalBoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q, GeoLib::Point const& r);
+ /// Bounding sphere using four points
+ MinimalBoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q, GeoLib::Point const& r, GeoLib::Point const& s);
+ /// Bounding sphere of n points
+ MinimalBoundingSphere(std::vector<GeoLib::Point*> const& points);
+ ~MinimalBoundingSphere() {}
+
+ /// Returns the center point of the sphere
+ GeoLib::Point getCenter() const { return GeoLib::Point(_center.getCoords()); }
+
+ /// Returns the radius of the sphere
+ double getRadius() const {return _radius; }
+
+ /// Returns the squared euclidean distance of a point from the sphere (for points within the sphere distance is negative)
+ double pointDistanceSquared(GeoLib::Point const& pnt) const;
+
+ /// Creates n_points random points located on the surface of the bounding sphere (useful for visualisation)
+ std::vector<GeoLib::Point*>* getRandomSpherePoints(std::size_t n_points) const;
+
+private:
+ /// Constructor using no points
+ MinimalBoundingSphere();
+
+ /**
+ * Recursive method for calculating a minimal bounding sphere for an arbitrary number of points.
+ * Note: This method will change the order of elements in the vector sphere_points.
+ * \param sphere_points The vector of points for which the smallest enclosing sphere is calculated
+ * \param start_idx Start index of the vector subrange analysed in the current recursive step
+ * \param length Length of the vector subrange analysed in the current recursive step
+ * \param n_boundary_points Number of found boundary points in the current recursive step
+ *
+ * Algorithm based the following two papers:
+ * Emo Welzl: Smallest enclosing disks (balls and ellipsoids). New Results and New Trends in Computer Science, pp. 359--370, 1991
+ * Bernd Gaertner: Fast and Robust Smallest Enclosing Balls. ESA99, pp. 325--338, 1999.
+ * Code based on "Smallest Enclosing Spheres" implementation by Nicolas Capens on flipcode's Developer Toolbox (www.flipcode.com)
+ */
+ static MinimalBoundingSphere recurseCalculation(std::vector<GeoLib::Point*> sphere_points, std::size_t start_idx, std::size_t length, std::size_t n_boundary_points);
+
+ double _radius;
+ MathLib::Vector3 _center;
+};
+
+} // namespace
+
+#endif /* MINIMALBOUNDINGSPHERE_H_ */
diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index a086660e996d0f8a01256e1c90eb0ba8a49489c7..fe2bbc9673ee0255eda98d9b30024b180370f637 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -62,11 +62,6 @@ double getAngle (const double p0[3], const double p1[3], const double p2[3])
return acos (scalarProduct<double,3> (v0,v1) / (sqrt(scalarProduct<double,3>(v0,v0)) * sqrt(scalarProduct<double,3>(v1,v1))));
}
-double calcTetrahedronVolume(const double* x1, const double* x2, const double* x3, const double* x4)
-{
- return fabs((x1[0] - x4[0]) * ((x2[1] - x4[1]) * (x3[2] - x4[2]) - (x2[2] - x4[2]) * (x3[1] - x4[1]))
- - (x1[1] - x4[1]) * ((x2[0] - x4[0]) * (x3[2] - x4[2]) - (x2[2] - x4[2]) * (x3[0] - x4[0]))
- + (x1[2] - x4[2]) * ((x2[0] - x4[0]) * (x3[1] - x4[1]) - (x2[1] - x4[1]) * (x3[0] - x4[0]))) / 6.0;
-}
+
} // namespace
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index 8c78efca61690d9c1f4cc4dca30118855648d0af..47357fc6368af9b61345ed03bbd8372ae00cf9de 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -148,11 +148,6 @@ float normalize(float min, float max, float val);
*/
double getAngle (const double p0[3], const double p1[3], const double p2[3]);
-/**
- * Calculates the volume of a tetrahedron.
- * The formula is V=1/6*|a(b x c)| with a=x1->x2, b=x1->x3 and c=x1->x4.
- */
-double calcTetrahedronVolume(const double* x1, const double* x2, const double* x3, const double* x4);
/**
* simple power function that takes as a second argument an integer instead of a float
diff --git a/MeshLib/Elements/TemplateHex-impl.h b/MeshLib/Elements/TemplateHex-impl.h
index a88e1f9d504c25a7b2660aceb9c5306ee6e9bea2..45e123cbd9c0e5004337beeb2bfe463f77b4e054 100644
--- a/MeshLib/Elements/TemplateHex-impl.h
+++ b/MeshLib/Elements/TemplateHex-impl.h
@@ -18,7 +18,7 @@
#include "Quad.h"
#include "Prism.h"
-#include "MathTools.h"
+#include "AnalyticalGeometry.h"
namespace MeshLib {
@@ -105,12 +105,12 @@ TemplateHex<NNODES,CELLHEXTYPE>::~TemplateHex()
template <unsigned NNODES, CellType CELLHEXTYPE>
double TemplateHex<NNODES,CELLHEXTYPE>::computeVolume()
{
- return MathLib::calcTetrahedronVolume(_nodes[4]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[0]->getCoords())
- + MathLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[3]->getCoords(), _nodes[1]->getCoords(), _nodes[0]->getCoords())
- + MathLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[7]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords())
- + MathLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[7]->getCoords(), _nodes[6]->getCoords(), _nodes[2]->getCoords())
- + MathLib::calcTetrahedronVolume(_nodes[1]->getCoords(), _nodes[3]->getCoords(), _nodes[5]->getCoords(), _nodes[2]->getCoords())
- + MathLib::calcTetrahedronVolume(_nodes[3]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[2]->getCoords());
+ return GeoLib::calcTetrahedronVolume(_nodes[4]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[0]->getCoords())
+ + GeoLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[3]->getCoords(), _nodes[1]->getCoords(), _nodes[0]->getCoords())
+ + GeoLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[7]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords())
+ + GeoLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[7]->getCoords(), _nodes[6]->getCoords(), _nodes[2]->getCoords())
+ + GeoLib::calcTetrahedronVolume(_nodes[1]->getCoords(), _nodes[3]->getCoords(), _nodes[5]->getCoords(), _nodes[2]->getCoords())
+ + GeoLib::calcTetrahedronVolume(_nodes[3]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[2]->getCoords());
}
template <unsigned NNODES, CellType CELLHEXTYPE>
diff --git a/MeshLib/Elements/TemplatePrism-impl.h b/MeshLib/Elements/TemplatePrism-impl.h
index 3029b24043788017370fa7653562f427b70547bf..a91d0a2e3c7e2459008a38401187aa1bc86d92b5 100644
--- a/MeshLib/Elements/TemplatePrism-impl.h
+++ b/MeshLib/Elements/TemplatePrism-impl.h
@@ -20,7 +20,7 @@
#include "Pyramid.h"
#include "Quad.h"
-#include "MathTools.h"
+#include "AnalyticalGeometry.h"
namespace MeshLib {
@@ -104,9 +104,9 @@ TemplatePrism<NNODES,CELLPRISMTYPE>::~TemplatePrism()
template <unsigned NNODES, CellType CELLPRISMTYPE>
double TemplatePrism<NNODES,CELLPRISMTYPE>::computeVolume()
{
- return MathLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
- + MathLib::calcTetrahedronVolume(_nodes[1]->getCoords(), _nodes[4]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
- + MathLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[4]->getCoords(), _nodes[5]->getCoords(), _nodes[3]->getCoords());
+ return GeoLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
+ + GeoLib::calcTetrahedronVolume(_nodes[1]->getCoords(), _nodes[4]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
+ + GeoLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[4]->getCoords(), _nodes[5]->getCoords(), _nodes[3]->getCoords());
}
template <unsigned NNODES, CellType CELLPRISMTYPE>
diff --git a/MeshLib/Elements/TemplatePyramid-impl.h b/MeshLib/Elements/TemplatePyramid-impl.h
index 4a93bab8bfc9a1e3404e64082ebbd0ef6bdc28d5..5bcdd4d0ac21eae981d7bdac57cc90459ce780e7 100644
--- a/MeshLib/Elements/TemplatePyramid-impl.h
+++ b/MeshLib/Elements/TemplatePyramid-impl.h
@@ -20,7 +20,7 @@
#include "Tet.h"
#include "Quad.h"
-#include "MathTools.h"
+#include "AnalyticalGeometry.h"
namespace MeshLib {
@@ -107,8 +107,8 @@ TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::~TemplatePyramid()
template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
double TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::computeVolume()
{
- return MathLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[4]->getCoords())
- + MathLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords(), _nodes[4]->getCoords());
+ return GeoLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[4]->getCoords())
+ + GeoLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords(), _nodes[4]->getCoords());
}
template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
diff --git a/MeshLib/Elements/TemplateTet-impl.h b/MeshLib/Elements/TemplateTet-impl.h
index 37a06afed399c9edffd91523583ba8a435aaa3e5..0f9b967dec8d05580ca9d7c3b891307a15328d7f 100644
--- a/MeshLib/Elements/TemplateTet-impl.h
+++ b/MeshLib/Elements/TemplateTet-impl.h
@@ -17,7 +17,7 @@
#include "Node.h"
#include "Tri.h"
-#include "MathTools.h"
+#include "AnalyticalGeometry.h"
namespace MeshLib {
@@ -99,7 +99,7 @@ TemplateTet<NNODES,CELLTETTYPE>::~TemplateTet()
template <unsigned NNODES, CellType CELLTETTYPE>
double TemplateTet<NNODES,CELLTETTYPE>::computeVolume()
{
- return MathLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords());
+ return GeoLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords());
}
template <unsigned NNODES, CellType CELLTETTYPE>
diff --git a/Tests/GeoLib/TestBoundingSphere.cpp b/Tests/GeoLib/TestBoundingSphere.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d2416bf7a9a4c0d986a1b46f483087b92bf2ee08
--- /dev/null
+++ b/Tests/GeoLib/TestBoundingSphere.cpp
@@ -0,0 +1,117 @@
+/**
+ * \file TestBoundingSphere.cpp
+ * \author Karsten Rink
+ * \date 2014-08-29
+ *
+ * \copyright
+ * Copyright (c) 2012-2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ *
+ */
+
+#include "gtest/gtest.h"
+
+#include "GeoLib/MinimalBoundingSphere.h"
+#include "GeoLib/Point.h"
+
+TEST(GeoLib, TestBoundingSphere)
+{
+ GeoLib::Point a(0, 0 ,0);
+ GeoLib::Point b(2, 0 ,0);
+ GeoLib::Point c(1, 0.1 ,0);
+ GeoLib::Point d(1, -0.1 ,0);
+ std::vector<GeoLib::Point*> pnts;
+ pnts.push_back(new GeoLib::Point(0, 0 , 0));
+ pnts.push_back(new GeoLib::Point(2, 0 , 0));
+ pnts.push_back(new GeoLib::Point(1, 0.1 , 0));
+ pnts.push_back(new GeoLib::Point(1, -0.1 , 0));
+
+ {
+ /**
+ * Four points located like this:
+ *
+ * *
+ * * *
+ * *
+ *
+ * Tests if a smaller number of points than available is used if the resulting sphere is smaller.
+ * Expected result is C=(1,0,0), r=1
+ */
+ GeoLib::MinimalBoundingSphere s(pnts);
+ GeoLib::Point center = s.getCenter();
+ ASSERT_NEAR(1.0, center[0], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.0, center[1], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.0, center[2], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(1.0, s.getRadius(), std::numeric_limits<double>::epsilon());
+ }
+
+ {
+ /**
+ * Four points located like this:
+ *
+ * *
+ * * *
+ *
+ *
+ * *
+ *
+ * The smallest sphere has a diameter that is larger than the distance between any two points.
+ * Expected result is C=(1,0.0246,-0.3446), r=1.058
+ */
+ (*pnts[2])[2] -= 1.4;
+ GeoLib::MinimalBoundingSphere s(pnts);
+ GeoLib::Point center = s.getCenter();
+ ASSERT_NEAR(1.0, center[0], 0.0001);
+ ASSERT_NEAR(0.0246, center[1], 0.0001);
+ ASSERT_NEAR(-0.3446, center[2], 0.0001);
+ ASSERT_NEAR(1.0580, s.getRadius(), 0.0001);
+ }
+
+ pnts[0] = new GeoLib::Point(0, 0, 0);
+ pnts[1] = new GeoLib::Point(1, 0, 0);
+ pnts[2] = new GeoLib::Point(1, 1, 0);
+ pnts[3] = new GeoLib::Point(0, 1, 0);
+ pnts.push_back(new GeoLib::Point(0, 0, 1));
+ pnts.push_back(new GeoLib::Point(1, 0, 1));
+ pnts.push_back(new GeoLib::Point(1, 1, 1));
+ pnts.push_back(new GeoLib::Point(0, 1, 0.9));
+
+ {
+ /**
+ * A "cube" where one node is pushed slightly towards the centre (and should be ignored).
+ * Expected result is C=(0.5,0.5,0.5), r=0.866
+ */
+ GeoLib::MinimalBoundingSphere s(pnts);
+ GeoLib::Point center = s.getCenter();
+ ASSERT_NEAR(0.5, center[0], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.5, center[1], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.5, center[2], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.8660, s.getRadius(), 0.0001);
+ }
+
+ /**
+ * A "cube" where one node is pulled away from the centre (making the resulting sphere larger).
+ * Expected result is C=(0.5,0.5,0.6), r=0.9273
+ */
+ (*pnts[7])[2] += 0.3;
+ GeoLib::MinimalBoundingSphere s(pnts);
+ {
+ GeoLib::Point center = s.getCenter();
+ ASSERT_NEAR(0.5, center[0], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.5, center[1], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.6, center[2], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.9273, s.getRadius(), 0.0001);
+ }
+
+ /// Calculates the bounding sphere of points on a bounding sphere
+ std::vector<GeoLib::Point*> *sphere_points (s.getRandomSpherePoints(1000));
+ GeoLib::MinimalBoundingSphere t(*sphere_points);
+ GeoLib::Point center = s.getCenter();
+ ASSERT_NEAR(0.5, center[0], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.5, center[1], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.6, center[2], std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.9273, s.getRadius(), 0.0001);
+}