diff --git a/GeoLib/BoundingSphere.cpp b/GeoLib/BoundingSphere.cpp new file mode 100644 index 0000000000000000000000000000000000000000..f99283ccac721bc2f535a24ce4b3686802569601 --- /dev/null +++ b/GeoLib/BoundingSphere.cpp @@ -0,0 +1,164 @@ +/** + * \file Calculation of a minimum bounding sphere for a vector of points + * \author Karsten Rink + * \date 2014-07-11 + * \brief Implementation of the BoundingSphere class. + * + * \copyright + * Copyright (c) 2013, 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 "BoundingSphere.h" + +// ThirdParty/logog +#include "logog/include/logog.hpp" + +#include "MathTools.h" + +namespace GeoLib { + +BoundingSphere::BoundingSphere() +: _center(0,0,0), _radius(-1) +{ +} + +BoundingSphere::BoundingSphere(const BoundingSphere &sphere) +: _center(sphere.getCenter()), _radius(sphere.getRadius()) +{ +} + +BoundingSphere::BoundingSphere(const GeoLib::Point &p) +: _center(p), _radius(std::numeric_limits<double>::epsilon()) +{ +} + +BoundingSphere::BoundingSphere(const GeoLib::Point &p, double radius) +: _center(p), _radius(radius) +{ +} + +BoundingSphere::BoundingSphere(const GeoLib::Point &p, const GeoLib::Point &q) +{ + const MathLib::Vector3 a(p, q); + const MathLib::Vector3 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) +{ + const MathLib::Vector3 a(p,r); + const MathLib::Vector3 b(p,q); + + const MathLib::Vector3 cross_ab(crossProduct(a,b)); + const double denom = 2.0 * scalarProduct(cross_ab,cross_ab); + const MathLib::Vector3 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; +} + +BoundingSphere::BoundingSphere(const GeoLib::Point &p, const GeoLib::Point &q, const GeoLib::Point &r, const GeoLib::Point &s) +{ + const MathLib::Vector3 a(p, q); + const MathLib::Vector3 b(p, r); + const MathLib::Vector3 c(p, 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]) + - 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) + + 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; +} + +BoundingSphere::BoundingSphere(const std::vector<GeoLib::Point*> &points) +: _center(0,0,0), _radius(-1) +{ + std::vector<GeoLib::Point*> sphere_points; + sphere_points.reserve(points.size()); + std::copy(points.cbegin(), points.cend(), std::back_inserter(sphere_points)); + + const BoundingSphere bounding_sphere = recurseCalculation(sphere_points, sphere_points.size(), 0); + this->_center = bounding_sphere.getCenter(); + this->_radius = bounding_sphere.getRadius(); +} + +BoundingSphere BoundingSphere::recurseCalculation(std::vector<GeoLib::Point*> &sphere_points, std::size_t idx, std::size_t boundary_points) +{ + BoundingSphere sphere; + switch(boundary_points) + { + case 0: + sphere = BoundingSphere(); + break; + case 1: + sphere = BoundingSphere(*sphere_points[0]); + break; + case 2: + sphere = BoundingSphere(*sphere_points[0], *sphere_points[1]); + break; + case 3: + sphere = BoundingSphere(*sphere_points[0], *sphere_points[1], *sphere_points[2]); + break; + case 4: + { + sphere = BoundingSphere(*sphere_points[0], *sphere_points[1], *sphere_points[2], *sphere_points[3]); + return sphere; + } + } + + for(std::size_t i=0; i<idx; ++i) + { + if(sphere.sqrPointDist(*sphere_points[i]) > 0) + { + for(std::size_t j=i; j>0; --j) + { + GeoLib::Point* tmp = sphere_points[j]; + sphere_points[j] = sphere_points[j-1]; + sphere_points[j - 1] = tmp; + } + + sphere = recurseCalculation(sphere_points, i, boundary_points+1); + } + } + return sphere; +} + +double BoundingSphere::sqrPointDist(const GeoLib::Point pnt) const +{ + return MathLib::sqrDist(_center.getCoords(), pnt.getCoords())-(_radius*_radius); +} + +std::vector<GeoLib::Point*>* BoundingSphere::getSpherePoints(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 fac (this->_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 new file mode 100644 index 0000000000000000000000000000000000000000..b39e371eb0892709382b2433215af40420763e54 --- /dev/null +++ b/GeoLib/BoundingSphere.h @@ -0,0 +1,73 @@ +/** + * \file Calculation of a minimum bounding sphere for a vector of points + * \author Karsten Rink + * \date 2014-07-11 + * \brief Definition of the BoundingSphere class. + * + * \copyright + * Copyright (c) 2013, 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 BOUNDINGSPHERE_H_ +#define BOUNDINGSPHERE_H_ + +#include <vector> + +#include "Vector3.h" +#include "Point.h" + +namespace GeoLib +{ + +class BoundingSphere +{ +public: + /// Constructor using no points + BoundingSphere(); + /// Copy constructor + BoundingSphere(const BoundingSphere &sphere); + /// Point-Sphere + BoundingSphere(const GeoLib::Point &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); + /// Bounding sphere of n points + BoundingSphere(const std::vector<GeoLib::Point*> &points); + ~BoundingSphere() {} + + /// 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 distance of a point from the sphere (for points within the sphere distance is negative) + double sqrPointDist(const GeoLib::Point pnt) const; + + std::vector<GeoLib::Point*>* getSpherePoints(std::size_t n_points) const; + +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. + * Code based on "Smallest Enclosing Spheres" by Nicolas Capens on flipcode's Developer Toolbox (www.flipcode.com) + */ + static BoundingSphere recurseCalculation(std::vector<GeoLib::Point*> &sphere_points, std::size_t idx, std::size_t boundary_points); + + double _radius; + MathLib::Vector3 _center; +}; + +} // namespace + +#endif /* BOUNDINGSPHERE_H_ */