Merge pull request #447 from rinkk/RadiusEdgeRatio

Minimal bounding spheres for points

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

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_ */

GeoLib/MinimalBoundingSphere.cpp

+/**
+ * \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;
+}
+
+}

GeoLib/MinimalBoundingSphere.h

+/**
+ * \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_ */

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

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

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>

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>

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>

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>

Tests/GeoLib/TestBoundingSphere.cpp

+/**
+ * \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);
+}