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;
18 BoundingSphere(): radius(0) { }
19 BoundingSphere(const LinAl::Vector<T, N> &, T);
21 template<typename Iter>
22 static BoundingSphere from_point_cloud(const Iter &, const Iter &);
24 const LinAl::Vector<T, N> get_center() const { return center; }
25 T get_radius() const { return radius; }
28 template<typename T, unsigned N>
29 BoundingSphere<T, N>::BoundingSphere(const LinAl::Vector<T, N> &c, T r):
34 template<typename T, unsigned N>
35 template<typename Iter>
36 BoundingSphere<T, N> BoundingSphere<T, N>::from_point_cloud(const Iter &begin, const Iter &end)
40 LinAl::Vector<T, N> p1;
42 for(Iter i=begin; i!=end; ++i)
44 T d = inner_product(*i, *i);
52 LinAl::Vector<T, N> p2;
54 for(Iter i=begin; i!=end; ++i)
56 LinAl::Vector<T, N> v = *i-p1;
57 T d = inner_product(v, v);
66 BoundingSphere<T, N> bsphere((p1+p2)/T(2), sqrt(sqdist));
67 for(Iter i=begin; i!=end; ++i)
69 LinAl::Vector<T, N> v = *i-bsphere.center;
70 T d = inner_product(v, v);
74 bsphere.center += v*(1-bsphere.radius/d);
75 bsphere.radius += d/2;
76 sqdist = bsphere.radius*bsphere.radius;
83 template<typename T, unsigned N>
84 BoundingSphere<T, N> operator|(const BoundingSphere<T, N> &bs1, const BoundingSphere<T, N> &bs2)
86 LinAl::Vector<T, N> v = bs2.get_center()-bs1.get_center();
88 if(d+bs2.get_radius()<bs1.get_radius())
90 else if(d+bs1.get_radius()<bs2.get_radius())
94 T r = (bs1.get_radius()+d+bs2.get_radius())/T(2);
95 if(bs1.get_radius()>bs2.get_radius())
96 return BoundingSphere<T, N>(bs1.get_center()+v*((r-bs1.get_radius())/d), r);
98 return BoundingSphere<T, N>(bs2.get_center()-v*((r-bs2.get_radius())/d), r);
102 } // namespace Geometry