X-Git-Url: http://git.tdb.fi/?a=blobdiff_plain;f=source%2Fgeometry%2Fboundingsphere.h;h=66c14d5c904776f68c7344e76e534e1638e4f933;hb=3c6fbd933b59fe247779146913fbca5b9656019f;hp=40826ca806915e529442c824bc290ce7140da819;hpb=69a0220ce9c71ebc92a4f9e77baf53780fee2206;p=libs%2Fmath.git

diff --git a/source/geometry/boundingsphere.h b/source/geometry/boundingsphere.h
index 40826ca..66c14d5 100644
--- 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())