Add a forgotten & to a getter

[libs/math.git] / source / geometry / boundingsphere.h

diff --git a/source/geometry/boundingsphere.h b/source/geometry/boundingsphere.h
index 40826ca806915e529442c824bc290ce7140da819..66c14d5c904776f68c7344e76e534e1638e4f933 100644 (file)
--- a/source/geometry/boundingsphere.h
+++ b/source/geometry/boundingsphere.h
@@ -13,22 +13,31 @@ class BoundingSphere
 private:
        LinAl::Vector<T, N> center;
        T radius;
+       bool empty;
 
 public:
-       BoundingSphere(): radius(0) { }
+       BoundingSphere();
        BoundingSphere(const LinAl::Vector<T, N> &, T);
 
        template<typename Iter>
        static BoundingSphere from_point_cloud(const Iter &, const Iter &);
 
-       const LinAl::Vector<T, N> get_center() const { return center; }
+       const LinAl::Vector<T, N> &get_center() const { return center; }
        T get_radius() const { return radius; }
+       bool is_empty() const { return empty; }
 };
 
+template<typename T, unsigned N>
+BoundingSphere<T, N>::BoundingSphere():
+       radius(0),
+       empty(true)
+{ }
+
 template<typename T, unsigned N>
 BoundingSphere<T, N>::BoundingSphere(const LinAl::Vector<T, N> &c, T r):
        center(c),
-       radius(r)
+       radius(r),
+       empty(false)
 { }
 
 template<typename T, unsigned N>
@@ -37,6 +46,9 @@ BoundingSphere<T, N> BoundingSphere<T, N>::from_point_cloud(const Iter &begin, c
 {
        using std::sqrt;
 
+       if(begin==end)
+               return BoundingSphere<T, N>();
+
        LinAl::Vector<T, N> p1;
        T sqdist = 0;
        for(Iter i=begin; i!=end; ++i)
@@ -71,8 +83,8 @@ BoundingSphere<T, N> BoundingSphere<T, N>::from_point_cloud(const Iter &begin, c
                if(d>sqdist)
                {
                        d = sqrt(d);
-                       bsphere.center += v*(1-bsphere.radius/d);
-                       bsphere.radius += d/2;
+                       bsphere.center += v*((T(1)-bsphere.radius/d)/T(2));
+                       bsphere.radius = (bsphere.radius+d)/T(2);
                        sqdist = bsphere.radius*bsphere.radius;
                }
        }
@@ -83,6 +95,11 @@ BoundingSphere<T, N> BoundingSphere<T, N>::from_point_cloud(const Iter &begin, c
 template<typename T, unsigned N>
 BoundingSphere<T, N> operator|(const BoundingSphere<T, N> &bs1, const BoundingSphere<T, N> &bs2)
 {
+       if(bs1.is_empty())
+               return bs2;
+       else if(bs2.is_empty())
+               return bs1;
+
        LinAl::Vector<T, N> v = bs2.get_center()-bs1.get_center();
        T d = v.norm();
        if(d+bs2.get_radius()<bs1.get_radius())