1 #ifndef MSP_GEOMETRY_BOUNDINGSPHERE_H_
2 #define MSP_GEOMETRY_BOUNDINGSPHERE_H_
5 #include <msp/linal/vector.h>
10 template<typename T, unsigned N>
14 LinAl::Vector<T, N> center;
20 BoundingSphere(const LinAl::Vector<T, N> &, T);
22 template<typename Iter>
23 static BoundingSphere from_point_cloud(const Iter &, const Iter &);
25 const LinAl::Vector<T, N> &get_center() const { return center; }
26 T get_radius() const { return radius; }
27 bool is_empty() const { return empty; }
30 template<typename T, unsigned N>
31 BoundingSphere<T, N>::BoundingSphere():
36 template<typename T, unsigned N>
37 BoundingSphere<T, N>::BoundingSphere(const LinAl::Vector<T, N> &c, T r):
43 template<typename T, unsigned N>
44 template<typename Iter>
45 BoundingSphere<T, N> BoundingSphere<T, N>::from_point_cloud(const Iter &begin, const Iter &end)
50 return BoundingSphere<T, N>();
52 LinAl::Vector<T, N> p1;
54 for(Iter i=begin; i!=end; ++i)
56 T d = inner_product(*i, *i);
64 LinAl::Vector<T, N> p2;
66 for(Iter i=begin; i!=end; ++i)
68 LinAl::Vector<T, N> v = *i-p1;
69 T d = inner_product(v, v);
78 BoundingSphere<T, N> bsphere((p1+p2)/T(2), sqrt(sqdist));
79 for(Iter i=begin; i!=end; ++i)
81 LinAl::Vector<T, N> v = *i-bsphere.center;
82 T d = inner_product(v, v);
86 bsphere.center += v*((T(1)-bsphere.radius/d)/T(2));
87 bsphere.radius = (bsphere.radius+d)/T(2);
88 sqdist = bsphere.radius*bsphere.radius;
95 template<typename T, unsigned N>
96 BoundingSphere<T, N> operator|(const BoundingSphere<T, N> &bs1, const BoundingSphere<T, N> &bs2)
100 else if(bs2.is_empty())
103 LinAl::Vector<T, N> v = bs2.get_center()-bs1.get_center();
105 if(d+bs2.get_radius()<bs1.get_radius())
107 else if(d+bs1.get_radius()<bs2.get_radius())
111 T r = (bs1.get_radius()+d+bs2.get_radius())/T(2);
112 if(bs1.get_radius()>bs2.get_radius())
113 return BoundingSphere<T, N>(bs1.get_center()+v*((r-bs1.get_radius())/d), r);
115 return BoundingSphere<T, N>(bs2.get_center()-v*((r-bs2.get_radius())/d), r);
119 } // namespace Geometry