From: Mikko Rasa Date: Wed, 11 Sep 2013 07:15:02 +0000 (+0300) Subject: Add a bounding sphere class X-Git-Url: http://git.tdb.fi/?a=commitdiff_plain;h=69a0220ce9c71ebc92a4f9e77baf53780fee2206;p=libs%2Fmath.git Add a bounding sphere class Shapes don't produce bounding spheres yet, but it can be used for other purposes. --- diff --git a/source/geometry/boundingsphere.h b/source/geometry/boundingsphere.h new file mode 100644 index 0000000..40826ca --- /dev/null +++ b/source/geometry/boundingsphere.h @@ -0,0 +1,105 @@ +#ifndef MSP_GEOMETRY_BOUNDINGSPHERE_H_ +#define MSP_GEOMETRY_BOUNDINGSPHERE_H_ + +#include +#include + +namespace Msp { +namespace Geometry { + +template +class BoundingSphere +{ +private: + LinAl::Vector center; + T radius; + +public: + BoundingSphere(): radius(0) { } + BoundingSphere(const LinAl::Vector &, T); + + template + static BoundingSphere from_point_cloud(const Iter &, const Iter &); + + const LinAl::Vector get_center() const { return center; } + T get_radius() const { return radius; } +}; + +template +BoundingSphere::BoundingSphere(const LinAl::Vector &c, T r): + center(c), + radius(r) +{ } + +template +template +BoundingSphere BoundingSphere::from_point_cloud(const Iter &begin, const Iter &end) +{ + using std::sqrt; + + LinAl::Vector p1; + T sqdist = 0; + for(Iter i=begin; i!=end; ++i) + { + T d = inner_product(*i, *i); + if(d>sqdist) + { + p1 = *i; + sqdist = d; + } + } + + LinAl::Vector p2; + sqdist = 0; + for(Iter i=begin; i!=end; ++i) + { + LinAl::Vector v = *i-p1; + T d = inner_product(v, v); + if(d>sqdist) + { + p2 = *i; + sqdist = d; + } + } + + sqdist /= 4; + BoundingSphere bsphere((p1+p2)/T(2), sqrt(sqdist)); + for(Iter i=begin; i!=end; ++i) + { + LinAl::Vector v = *i-bsphere.center; + T d = inner_product(v, v); + if(d>sqdist) + { + d = sqrt(d); + bsphere.center += v*(1-bsphere.radius/d); + bsphere.radius += d/2; + sqdist = bsphere.radius*bsphere.radius; + } + } + + return bsphere; +} + +template +BoundingSphere operator|(const BoundingSphere &bs1, const BoundingSphere &bs2) +{ + LinAl::Vector v = bs2.get_center()-bs1.get_center(); + T d = v.norm(); + if(d+bs2.get_radius()bs2.get_radius()) + return BoundingSphere(bs1.get_center()+v*((r-bs1.get_radius())/d), r); + else + return BoundingSphere(bs2.get_center()-v*((r-bs2.get_radius())/d), r); + } +} + +} // namespace Geometry +} // namespace Msp + +#endif