diff --git a/GeoLib/AABB.h b/GeoLib/AABB.h
index 36a100bebb8720169de45f2bc87d0a6c12ead264..b5ed9b59e3062e03effc4c728653178c7e8e7d17 100644
--- a/GeoLib/AABB.h
+++ b/GeoLib/AABB.h
@@ -15,11 +15,14 @@
 #ifndef AABB_H_
 #define AABB_H_
 
-#include <limits>
+#include <bitset>
+#include <cassert>
+#include <cmath>
 #include <cstddef>
 #include <cstdlib>
 #include <iterator>
-#include <cassert>
+#include <limits>
+#include <tuple>
 #include <vector>
 
 #include <logog/include/logog.hpp>
@@ -49,8 +52,9 @@ public:
 		assert(! ids.empty());
 		init(pnts[ids[0]]);
 		for (std::size_t i=1; i<ids.size(); ++i) {
-			update(*(pnts[ids[i]]));
+			updateWithoutEnlarge(*(pnts[ids[i]]));
 		}
+		enlarge();
 	}
 
 	/**
@@ -81,25 +85,41 @@ public:
 		init(*first);
 		InputIterator it(first);
 		while (it != last) {
-			update(*it);
+			updateWithoutEnlarge(*it);
 			it++;
 		}
+		enlarge();
 	}
 
-	bool update(PNT_TYPE const & pnt)
+	/// Checks if the bounding box has to be updated.
+	/// @return true if AABB is updated.
+	bool update(PNT_TYPE const & p)
 	{
-		bool updated(false);
+		// First component of the pair signals if the minimum point is changed
+		// Second component signals not only if the max point is changed.
+		// Furthermore it is signaled what coordinate (0,1,2) is changed.
+		std::pair<bool,std::bitset<3>> updated(0,0);
 		for (std::size_t k(0); k<3; k++) {
-			if (pnt[k] < _min_pnt[k]) {
-				_min_pnt[k] = pnt[k];
-				updated = true;
+			// if the minimum point is updated pair.first==true
+			if (p[k] < _min_pnt[k]) {
+				_min_pnt[k] = p[k];
+				updated.first = true;
 			}
-			if (_max_pnt[k] < pnt[k]) {
-				_max_pnt[k] = pnt[k];
-				updated = true;
+			// if the kth coordinate of the maximum point is updated
+			// pair.second[k]==true
+			if (p[k] >= _max_pnt[k]) {
+				_max_pnt[k] = p[k];
+				updated.second[k] = true;
 			}
 		}
-		return updated;
+
+		if (updated.second.any()) {
+			enlarge(updated.second);
+			return true;
+		} else if (updated.first) {
+			return true;
+		}
+		return false;
 	}
 
 	/**
@@ -150,6 +170,19 @@ protected:
 		std::numeric_limits<double>::lowest(),
 		std::numeric_limits<double>::lowest()}}};
 private:
+	/// Enlarge the bounding box the smallest possible amount (modifying the
+	/// unit in the last place). Only the coordinates of the maximum point are
+	/// changed such that the half-open property will be preserved.
+	void enlarge(std::bitset<3> to_update = 7)
+	{
+		for (std::size_t k=0; k<3; ++k) {
+			if (to_update[k]) {
+				_max_pnt[k] = std::nextafter(_max_pnt[k],
+					std::numeric_limits<double>::max());
+			}
+		}
+	}
+
 	void init(PNT_TYPE const & pnt)
 	{
 		_min_pnt[0] = _max_pnt[0] = pnt[0];
@@ -160,9 +193,32 @@ private:
 	{
 		init(*pnt);
 	}
+
+	/// Private method that is used internally to update the min and max point
+	/// of the bounding box using point \f$p\f$ without enlarging the bounding
+	/// box. Using this method the bounding box of the initial point set is
+	/// enlarged only once.
+	/// @param p point that will possibly change the bounding box points
+	void  updateWithoutEnlarge(PNT_TYPE const & p)
+	{
+		for (std::size_t k(0); k<3; k++) {
+			if (p[k] < _min_pnt[k]) {
+				_min_pnt[k] = p[k];
+			}
+			if (p[k] >= _max_pnt[k]) {
+				_max_pnt[k] = p[k];
+			}
+		}
+	}
+
+	void updateWithoutEnlarge(PNT_TYPE const * pnt)
+	{
+		updateWithoutEnlarge(*pnt);
+	}
+
 	void update(PNT_TYPE const * pnt)
 	{
-		update (*pnt);
+		update(*pnt);
 	}
 };
 } // end namespace