diff --git a/GeoLib/BoundingSphere.cpp b/GeoLib/BoundingSphere.cpp
index 3818624019235b0de755c7f05738f73e22387a27..5d81f3903deb7502ae6fbced2fdc8cb30ad941e8 100644
--- a/GeoLib/BoundingSphere.cpp
+++ b/GeoLib/BoundingSphere.cpp
@@ -14,8 +14,7 @@
 
 #include "BoundingSphere.h"
 
-// ThirdParty/logog
-#include "logog/include/logog.hpp"
+#include <ctime>
 
 #include "MathTools.h"
 #include "AnalyticalGeometry.h"
@@ -27,45 +26,45 @@ BoundingSphere::BoundingSphere()
 {	
 }
 
-BoundingSphere::BoundingSphere(const BoundingSphere &sphere)
+BoundingSphere::BoundingSphere(BoundingSphere const& sphere)
 : _center(sphere.getCenter()), _radius(sphere.getRadius())
 {
 }
 
-BoundingSphere::BoundingSphere(const GeoLib::Point &p)
+BoundingSphere::BoundingSphere(GeoLib::Point const& p)
 : _center(p), _radius(std::numeric_limits<double>::epsilon())
 {
 }
 
-BoundingSphere::BoundingSphere(const GeoLib::Point &p, double radius)
+BoundingSphere::BoundingSphere(GeoLib::Point const& p, double radius)
 : _center(p), _radius(radius)
 {
 }
 
-BoundingSphere::BoundingSphere(const GeoLib::Point &p, const GeoLib::Point &q)
+BoundingSphere::BoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q)
 : _center(p), _radius(std::numeric_limits<double>::epsilon())
 {
-	const MathLib::Vector3 a(p, q);
+	MathLib::Vector3 const a(p, q);
 
     if (a.getLength() > 0)
     {
-	    const MathLib::Vector3 o(0.5*a);
+	    MathLib::Vector3 const o(0.5*a);
 	    _radius = o.getLength() + std::numeric_limits<double>::epsilon();
 	    _center = MathLib::Vector3(p) + o;
     }
 }
 
-BoundingSphere::BoundingSphere(const GeoLib::Point &p, const GeoLib::Point &q,  const GeoLib::Point &r)
+BoundingSphere::BoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q, GeoLib::Point const& r)
 {
-	const MathLib::Vector3 a(p,r);
-	const MathLib::Vector3 b(p,q);
+	MathLib::Vector3 const a(p,r);
+	MathLib::Vector3 const b(p,q);
 
-	const MathLib::Vector3 cross_ab(crossProduct(a,b));
+	MathLib::Vector3 const cross_ab(crossProduct(a,b));
 
     if (cross_ab.getLength() > 0)
     {
-	    const double denom = 2.0 * scalarProduct(cross_ab,cross_ab);
-	    const MathLib::Vector3 o = (scalarProduct(b,b) * crossProduct(cross_ab, a) 
+	    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();
@@ -83,32 +82,32 @@ BoundingSphere::BoundingSphere(const GeoLib::Point &p, const GeoLib::Point &q,
     }
 }
 
-BoundingSphere::BoundingSphere(const GeoLib::Point &p, const GeoLib::Point &q, const GeoLib::Point &r, const GeoLib::Point &s)
+BoundingSphere::BoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q, GeoLib::Point const& r, GeoLib::Point const& s)
 {
-	const MathLib::Vector3 a(p, q);
-	const MathLib::Vector3 b(p, r);
-	const MathLib::Vector3 c(p, 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))
     {
     	// det of matrix [a^T, b^T, c^T]^T
-	    const double denom = 2.0 * (a[0] * (b[1] * c[2] - c[1] * b[2])
+	    double const denom = 2.0 * (a[0] * (b[1] * c[2] - c[1] * b[2])
 							      - b[0] * (a[1] * c[2] - c[1] * a[2])
 							      + c[0] * (a[1] * b[2] - b[1] * a[2]));
-	    const MathLib::Vector3 o = (scalarProduct(c,c) * crossProduct(a,b) 
+	    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);
+							      * (1.0 / denom);
 
 	    _radius = o.getLength() + std::numeric_limits<double>::epsilon();
 	    _center = MathLib::Vector3(p) + o;
     }
     else
     {
-        BoundingSphere pqr(p, q , r);
-        BoundingSphere pqs(p, q , s);
-        BoundingSphere prs(p, r , s);
-        BoundingSphere qrs(q, r , s);
+        BoundingSphere const pqr(p, q , r);
+        BoundingSphere const pqs(p, q , s);
+        BoundingSphere const prs(p, r , s);
+        BoundingSphere const qrs(q, r , s);
         _radius = pqr.getRadius();
         _center = pqr.getCenter();
         if (_radius < pqs.getRadius())
@@ -129,22 +128,18 @@ BoundingSphere::BoundingSphere(const GeoLib::Point &p, const GeoLib::Point &q, c
     }
 }
 
-BoundingSphere::BoundingSphere(const std::vector<GeoLib::Point*> &points)
+BoundingSphere::BoundingSphere(std::vector<GeoLib::Point*> const& points)
 : _center(0,0,0), _radius(-1)
 {
-	const std::size_t n_points (points.size());
-	GeoLib::Point **sphere_points = new GeoLib::Point*[n_points];
-	for(unsigned int i = 0; i < n_points; i++)
-		sphere_points[i] = points[i];
-
-	const BoundingSphere bounding_sphere = recurseCalculation(points, 0, points.size(), 0);
-	delete[] sphere_points;
-	
-	this->_center = bounding_sphere.getCenter();
-	this->_radius = bounding_sphere.getRadius();
+	std::size_t const n_points (points.size());
+	std::vector<GeoLib::Point*> sphere_points(points);
+ 
+	BoundingSphere const bounding_sphere = recurseCalculation(sphere_points, 0, sphere_points.size(), 0);
+	_center = bounding_sphere.getCenter();
+	_radius = bounding_sphere.getRadius();
 }
 
-BoundingSphere BoundingSphere::recurseCalculation(std::vector<GeoLib::Point*> sphere_points, std::size_t current_index, std::size_t n_points, std::size_t n_boundary_points)
+BoundingSphere BoundingSphere::recurseCalculation(std::vector<GeoLib::Point*> sphere_points, std::size_t start_idx, std::size_t length, std::size_t n_boundary_points)
 {
 	BoundingSphere sphere;
 	switch(n_boundary_points)
@@ -153,35 +148,37 @@ BoundingSphere BoundingSphere::recurseCalculation(std::vector<GeoLib::Point*> sp
 		sphere = BoundingSphere();
 		break;
 	case 1:
-		sphere = BoundingSphere(*sphere_points[current_index-1]);
+		sphere = BoundingSphere(*sphere_points[start_idx-1]);
 		break;
 	case 2:
-		sphere = BoundingSphere(*sphere_points[current_index-1], *sphere_points[current_index-2]);
+		sphere = BoundingSphere(*sphere_points[start_idx-1], *sphere_points[start_idx-2]);
 		break;
 	case 3:
-		sphere = BoundingSphere(*sphere_points[current_index-1], *sphere_points[current_index-2], *sphere_points[current_index-3]);
+        sphere = BoundingSphere(*sphere_points[start_idx-1], *sphere_points[start_idx-2], *sphere_points[start_idx-3]);
 		break;
 	case 4:
-		sphere = BoundingSphere(*sphere_points[current_index-1], *sphere_points[current_index-2], *sphere_points[current_index-3], *sphere_points[current_index-4]);
+        sphere = BoundingSphere(*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=current_index; i<n_points; ++i)
+	for(std::size_t i=0; i<length; ++i)
 	{
         // current point is located outside of sphere
-		if(sphere.sqrPointDist(*sphere_points[i]) > 0)
-		{
-			GeoLib::Point* tmp = sphere_points[i];
-            std::copy(sphere_points.begin(), sphere_points.begin() + i, sphere_points.begin() + 1);
-            sphere_points[0] = tmp;
-
-    		sphere = recurseCalculation(sphere_points, current_index+1, i, n_boundary_points+1);
+        if (sphere.sqrPointDist(*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 BoundingSphere::sqrPointDist(const GeoLib::Point pnt) const
+double BoundingSphere::sqrPointDist(GeoLib::Point const& pnt) const
 {
 	return MathLib::sqrDist(_center.getCoords(), pnt.getCoords())-(_radius*_radius);
 }
@@ -201,7 +198,7 @@ std::vector<GeoLib::Point*>* BoundingSphere::getRandomSpherePoints(std::size_t n
 			vec[i] = (double)rand()-(RAND_MAX/2.0);
 			sum+=(vec[i]*vec[i]);
 		}
-		double fac (this->_radius/sqrt(sum));
+		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/BoundingSphere.h b/GeoLib/BoundingSphere.h
index 711986b8a205e46dbb55220064fdd1b8faaf32fe..655bf3d307a450e59fa0308e5a364590152a33a6 100644
--- a/GeoLib/BoundingSphere.h
+++ b/GeoLib/BoundingSphere.h
@@ -30,19 +30,19 @@ public:
 	/// Constructor using no points
 	BoundingSphere();
 	/// Copy constructor
-	BoundingSphere(const BoundingSphere &sphere);
+	BoundingSphere(BoundingSphere const& sphere);
 	/// Point-Sphere
-	BoundingSphere(const GeoLib::Point &p);
+	BoundingSphere(GeoLib::Point const& p);
 	/// Constructor using center and radius
-	BoundingSphere(const GeoLib::Point &p, double radius);
-	/// Sphere through two points
-	BoundingSphere(const GeoLib::Point &p, const GeoLib::Point &q);
-	/// Sphere through three points
-	BoundingSphere(const GeoLib::Point &p, const GeoLib::Point &q, const GeoLib::Point &);
-	/// Sphere through four points
-	BoundingSphere(const GeoLib::Point &p, const GeoLib::Point &q, const GeoLib::Point &R, const Point &S);
+	BoundingSphere(GeoLib::Point const& p, double radius);
+	/// Bounding sphere using two points
+	BoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q);
+	/// Bounding sphere using three points
+	BoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q, GeoLib::Point const& r);
+	/// Bounding sphere using four points
+	BoundingSphere(GeoLib::Point const& p, GeoLib::Point const& q, GeoLib::Point const& r, GeoLib::Point const& s);
 	/// Bounding sphere of n points
-	BoundingSphere(const std::vector<GeoLib::Point*> &points);
+	BoundingSphere(std::vector<GeoLib::Point*> const& points);
 	~BoundingSphere() {}
 
 	/// Returns the center point of the sphere
@@ -52,7 +52,7 @@ public:
 	double getRadius() const {return _radius; }
 
 	/// Returns the squared distance of a point from the sphere (for points within the sphere distance is negative)
-	double sqrPointDist(const GeoLib::Point pnt) const;
+	double sqrPointDist(GeoLib::Point const& pnt) const;
 
 	/// Creates n_points random points located on the surface of the sphere (useful for visualisation)
 	std::vector<GeoLib::Point*>* getRandomSpherePoints(std::size_t n_points) const;
@@ -60,10 +60,18 @@ public:
 private:
 	/**
 	 * Recursive method for calculating a minimal bounding sphere for an arbitrary number of points.
-	 * Algorithm based on Bernd Gärtner: Fast and Robust Smallest Enclosing Balls. ESA99, pages 325-338, 1999.
+     * 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 Gärtner: 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 BoundingSphere recurseCalculation(std::vector<GeoLib::Point*> sphere_points, std::size_t current_index, std::size_t n_points, std::size_t n_boundary_points);
+	static BoundingSphere 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;