diff --git a/GeoLib/AABB.h b/GeoLib/AABB.h
index a9735afb9fc2c7ba530674a743afaedb07221f52..6b9563be5187922c1788c49b68d9fff4f7d3037e 100644
--- a/GeoLib/AABB.h
+++ b/GeoLib/AABB.h
@@ -13,9 +13,13 @@
 #ifndef AABB_H_
 #define AABB_H_
 
-#include "Point.h"
 #include <limits>
 #include <cstddef>
+#include <iterator>
+#include <cassert>
+#include <vector>
+
+#include "Point.h"
 
 namespace GeoLib
 {
@@ -33,10 +37,15 @@ public:
 	/**
 	 * construction of object, initialization the axis aligned bounding box
 	 * */
-	AABB() :
+	AABB(std::vector<PNT_TYPE*> const& pnts, std::vector<std::size_t> const& ids) :
 		_min_pnt(std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()),
 		_max_pnt(std::numeric_limits<double>::min(), std::numeric_limits<double>::min(), std::numeric_limits<double>::min())
-	{}
+	{
+		assert(! ids.empty());
+		for (auto i : ids) {
+			update(*(pnts[i]));
+		}
+	}
 
 	/**
 	 * copy constructor.
@@ -60,6 +69,8 @@ public:
 		_min_pnt(std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()),
 		_max_pnt(std::numeric_limits<double>::min(), std::numeric_limits<double>::min(), std::numeric_limits<double>::min())
 	{
+		assert(std::distance(first,last) > 0);
+		init(*first);
 		InputIterator it(first);
 		while (it != last) {
 			update(*it);
@@ -118,6 +129,16 @@ protected:
 	PNT_TYPE _min_pnt;
 	PNT_TYPE _max_pnt;
 private:
+	void init(PNT_TYPE const & pnt)
+	{
+		_min_pnt[0] = _max_pnt[0] = pnt[0];
+		_min_pnt[1] = _max_pnt[1] = pnt[1];
+		_min_pnt[2] = _max_pnt[2] = pnt[2];
+	}
+	void init(PNT_TYPE const * pnt)
+	{
+		init(*pnt);
+	}
 	void update(PNT_TYPE const * pnt)
 	{
 		update (*pnt);
diff --git a/GeoLib/PointVec.cpp b/GeoLib/PointVec.cpp
index 7937e6c4f420b9bb22827e06a83c58cd988e2960..a9f1e8c03f092d4791adbd07c075b8b491764f0e 100644
--- a/GeoLib/PointVec.cpp
+++ b/GeoLib/PointVec.cpp
@@ -26,12 +26,11 @@ namespace GeoLib
 PointVec::PointVec (const std::string& name, std::vector<Point*>* points,
                     std::map<std::string, size_t>* name_id_map, PointType type, double rel_eps) :
                     	TemplateVec<Point> (name, points, name_id_map),
-	_type(type), _sqr_shortest_dist (std::numeric_limits<double>::max())
+	_type(type), _sqr_shortest_dist (std::numeric_limits<double>::max()), _aabb(points->begin(), points->end())
 {
 	assert (_data_vec);
 	size_t number_of_all_input_pnts (_data_vec->size());
 
-	calculateAABB ();
 	rel_eps *= sqrt(MathLib::sqrDist (&(_aabb.getMinPoint()),&(_aabb.getMaxPoint())));
 	makePntsUnique (_data_vec, _pnt_id_map, rel_eps);
 
@@ -222,13 +221,6 @@ void PointVec::calculateShortestDistance ()
 	_sqr_shortest_dist = MathLib::sqrDist ((*_data_vec)[i], (*_data_vec)[j]);
 }
 
-void PointVec::calculateAABB ()
-{
-	const size_t n_pnts (_data_vec->size());
-	for (size_t i(0); i < n_pnts; i++)
-		_aabb.update (*(*_data_vec)[i]);
-}
-
 std::vector<GeoLib::Point*>* PointVec::getSubset(const std::vector<size_t> &subset)
 {
 	std::vector<GeoLib::Point*> *new_points (new std::vector<GeoLib::Point*>(subset.size()));
diff --git a/GeoLib/PointVec.h b/GeoLib/PointVec.h
index 76447841024003e286704717602b2c743c220c15..16eb08f6b7bdc1534c679683ae6f0a32a3aa41de 100644
--- a/GeoLib/PointVec.h
+++ b/GeoLib/PointVec.h
@@ -139,7 +139,6 @@ private:
 	 */
 	double _sqr_shortest_dist;
 
-	void calculateAABB ();
 	AABB<GeoLib::Point> _aabb;
 };
 } // end namespace
diff --git a/GeoLib/Polygon.cpp b/GeoLib/Polygon.cpp
index 835a1f4b2cd5868ab9752dcefd5b604c30641983..5f4d4bc833c3fd4db5e45689e4d23f9f3e634a02 100644
--- a/GeoLib/Polygon.cpp
+++ b/GeoLib/Polygon.cpp
@@ -29,16 +29,12 @@
 namespace GeoLib
 {
 Polygon::Polygon(const Polyline &ply, bool init) :
-	Polyline (ply)
+	Polyline(ply), _aabb(ply.getPointsVec(), ply._ply_pnt_ids)
 {
 	if (init)
 		initialise ();
 }
 
-Polygon::Polygon (const std::vector<Point*>& pnt_vec) :
-	Polyline (pnt_vec)
-{}
-
 Polygon::~Polygon()
 {
 	// remove polygons from list
@@ -52,7 +48,6 @@ Polygon::~Polygon()
 bool Polygon::initialise ()
 {
 	if (this->isClosed()) {
-		calculateAABB();
 		ensureCWOrientation();
 		return true;
 	} else {
@@ -234,13 +229,6 @@ EdgeType::value Polygon::getEdgeType (size_t k, GeoLib::Point const & pnt) const
 	}
 }
 
-void Polygon::calculateAABB ()
-{
-	size_t n_nodes (getNumberOfPoints());
-	for (size_t k(0); k < n_nodes; k++)
-		_aabb.update ((*(getPoint(k))));
-}
-
 void Polygon::ensureCWOrientation ()
 {
 	// *** pre processing: rotate points to xy-plan
@@ -398,14 +386,14 @@ void Polygon::splitPolygonAtPoint (std::list<GeoLib::Polygon*>::iterator polygon
 				std::swap (idx0, idx1);
 
 			// create two closed polylines
-			GeoLib::Polygon* polygon0 (new GeoLib::Polygon((*polygon_it)->getPointsVec()));
+			GeoLib::Polygon* polygon0 (new GeoLib::Polygon(*(*polygon_it)));
 			for (size_t k(0); k <= idx0; k++)
 				polygon0->addPoint ((*polygon_it)->getPointID (k));
 			for (size_t k(idx1 + 1); k < (*polygon_it)->getNumberOfPoints(); k++)
 				polygon0->addPoint ((*polygon_it)->getPointID (k));
 			polygon0->initialise();
 
-			GeoLib::Polygon* polygon1 (new GeoLib::Polygon((*polygon_it)->getPointsVec()));
+			GeoLib::Polygon* polygon1 (new GeoLib::Polygon(*(*polygon_it)));
 			for (size_t k(idx0); k <= idx1; k++)
 				polygon1->addPoint ((*polygon_it)->getPointID (k));
 			polygon1->initialise();
diff --git a/GeoLib/Polygon.h b/GeoLib/Polygon.h
index 3afde0f02e441fd4f020c75edd2768ac6f493ccb..87ba4e24e503b5220a5eba5ebcf4a7564c62fd40 100644
--- a/GeoLib/Polygon.h
+++ b/GeoLib/Polygon.h
@@ -53,8 +53,6 @@ public:
 	 */
 	Polygon(const Polyline &ply, bool init = true);
 
-	Polygon (const std::vector<Point*>& pnt_vec);
-
 	virtual ~Polygon();
 
 	/**
@@ -121,7 +119,6 @@ private:
 	 */
 	EdgeType::value getEdgeType (std::size_t k, GeoLib::Point const & pnt) const;
 
-	void calculateAABB ();
 	void ensureCWOrientation ();
 
 	void splitPolygonAtIntersection (std::list<Polygon*>::iterator polygon_it);
diff --git a/GeoLib/Polyline.h b/GeoLib/Polyline.h
index 1d3748c2919b850963c7577e79b5aca2e4abc47a..081c3666c087a9c711fbd0081a0d90f26b3ddf9e 100644
--- a/GeoLib/Polyline.h
+++ b/GeoLib/Polyline.h
@@ -50,6 +50,7 @@ public:
 class Polyline : public GeoObject
 {
 public:
+	friend class Polygon;
 	/** constructor
 	 * \param pnt_vec a reference to the point vector
 	 */
@@ -58,7 +59,7 @@ public:
 	 * Copy constructor
 	 * @param ply Polyline
 	 */
-	Polyline (const Polyline& ply);
+	Polyline(const Polyline& ply);
 
 	virtual ~Polyline() {}
 
diff --git a/GeoLib/Surface.cpp b/GeoLib/Surface.cpp
index 729923e852f77892a13066e0aa18efbe5c273e20..9675acc7758ffe51351d50f319ecb11e17edb3b4 100644
--- a/GeoLib/Surface.cpp
+++ b/GeoLib/Surface.cpp
@@ -24,16 +24,16 @@
 namespace GeoLib {
 
 Surface::Surface (const std::vector<Point*> &pnt_vec) :
-	GeoObject(), _sfc_pnts(pnt_vec), _bv()
+	GeoObject(), _sfc_pnts(pnt_vec), _bv(nullptr)
 {}
 
 Surface::~Surface ()
 {
-	for (size_t k(0); k<_sfc_triangles.size(); k++)
+	for (std::size_t k(0); k<_sfc_triangles.size(); k++)
 		delete _sfc_triangles[k];
 }
 
-void Surface::addTriangle (size_t pnt_a, size_t pnt_b, size_t pnt_c)
+void Surface::addTriangle (std::size_t pnt_a, std::size_t pnt_b, std::size_t pnt_c)
 {
 	assert (pnt_a < _sfc_pnts.size() && pnt_b < _sfc_pnts.size() && pnt_c < _sfc_pnts.size());
 
@@ -41,9 +41,13 @@ void Surface::addTriangle (size_t pnt_a, size_t pnt_b, size_t pnt_c)
 	if (pnt_a == pnt_b || pnt_a == pnt_c || pnt_b == pnt_c) return;
 
 	_sfc_triangles.push_back (new Triangle(_sfc_pnts, pnt_a, pnt_b, pnt_c));
-	_bv.update (*_sfc_pnts[pnt_a]);
-	_bv.update (*_sfc_pnts[pnt_b]);
-	_bv.update (*_sfc_pnts[pnt_c]);
+	if (!_bv) {
+		_bv = new AABB<Point>(_sfc_pnts, {pnt_a, pnt_b, pnt_c});
+	} else {
+		_bv->update (*_sfc_pnts[pnt_a]);
+		_bv->update (*_sfc_pnts[pnt_b]);
+		_bv->update (*_sfc_pnts[pnt_c]);
+	}
 }
 
 Surface* Surface::createSurface(const Polyline &ply)
@@ -86,12 +90,12 @@ Surface* Surface::createSurface(const Polyline &ply)
 
 }
 
-size_t Surface::getNTriangles () const
+std::size_t Surface::getNTriangles () const
 {
 	return _sfc_triangles.size();
 }
 
-const Triangle* Surface::operator[] (size_t i) const
+const Triangle* Surface::operator[] (std::size_t i) const
 {
 	assert (i < _sfc_triangles.size());
 	return _sfc_triangles[i];
@@ -99,13 +103,13 @@ const Triangle* Surface::operator[] (size_t i) const
 
 bool Surface::isPntInBV (const double *pnt) const
 {
-	return _bv.containsPoint (pnt);
+	return _bv->containsPoint (pnt);
 }
 
 bool Surface::isPntInSfc (const double *pnt) const
 {
 	bool nfound (true);
-	for (size_t k(0); k<_sfc_triangles.size() && nfound; k++) {
+	for (std::size_t k(0); k<_sfc_triangles.size() && nfound; k++) {
 		if (_sfc_triangles[k]->containsPoint (pnt)) {
 			nfound = false;
 		}
diff --git a/GeoLib/Surface.h b/GeoLib/Surface.h
index 262ba3d905cd6a0e49cb79feea0b93a15661920f..ac839a08f712179afbdc223bcbdb39579346be79 100644
--- a/GeoLib/Surface.h
+++ b/GeoLib/Surface.h
@@ -71,7 +71,7 @@ public:
 	 * method allows access to the internal axis aligned bounding box
 	 * @return axis aligned bounding box
 	 */
-	AABB<GeoLib::Point> const & getAABB () const { return _bv; }
+	AABB<GeoLib::Point> const& getAABB () const { return *_bv; }
 
 protected:
 	/** a vector of pointers to Points */
@@ -79,7 +79,7 @@ protected:
 	/** position of pointers to the geometric points */
 	std::vector<Triangle*> _sfc_triangles;
 	/** bounding volume is an axis aligned bounding box */
-	AABB<GeoLib::Point> _bv;
+	AABB<GeoLib::Point> *_bv;
 };
 
 }
diff --git a/Tests/GeoLib/TestAABB.cpp b/Tests/GeoLib/TestAABB.cpp
index dbae3c570fed931a78b9ae36acbcefc11ef2217d..85794d725c235824e1674db040e0f40fa323dbe3 100644
--- a/Tests/GeoLib/TestAABB.cpp
+++ b/Tests/GeoLib/TestAABB.cpp
@@ -21,14 +21,14 @@
 #include "Point.h"
 #include "AABB.h"
 
-TEST(GeoLibAABB, RandomNumberOfPointsPointersToRandomPointsInAList)
+TEST(GeoLibAABB, RandomNumberOfPointersToRandomPoints)
 {
 	/* initialize random seed: */
 	 srand ( time(NULL) );
 	 int n (rand() % 100000);
 	 int box_size (rand());
 	 int half_box_size(box_size/2);
-	 int minus_half_box_size(-1.0 * half_box_size);
+	 int minus_half_box_size(-half_box_size);
 
 	 // fill list with points
 	 std::list<GeoLib::Point*> pnts_list;
@@ -36,8 +36,6 @@ TEST(GeoLibAABB, RandomNumberOfPointsPointersToRandomPointsInAList)
 		 pnts_list.push_back(new GeoLib::Point(rand() % box_size - half_box_size, rand() % box_size - half_box_size, rand() % box_size - half_box_size));
 	 }
 
-	 std::cout << "testing with " << n << " points" << std::endl;
-
 	 // construct from list points a axis algined bounding box
 	 GeoLib::AABB<GeoLib::Point> aabb(pnts_list.begin(), pnts_list.end());
 
@@ -63,7 +61,7 @@ TEST(GeoLibAABB, RandomNumberOfPointsRandomPointInAList)
 	 int n (rand() % 1000000);
 	 int box_size (rand());
 	 int half_box_size(box_size/2);
-	 int minus_half_box_size(-1.0 * half_box_size);
+	 int minus_half_box_size(-half_box_size);
 
 	 // fill list with points
 	 std::list<GeoLib::Point> pnts_list;
@@ -71,8 +69,6 @@ TEST(GeoLibAABB, RandomNumberOfPointsRandomPointInAList)
 		 pnts_list.push_back(GeoLib::Point(rand() % box_size - half_box_size, rand() % box_size - half_box_size, rand() % box_size - half_box_size));
 	 }
 
-	 std::cout << "testing with " << n << " points" << std::endl;
-
 	 // construct from list points a axis algined bounding box
 	 GeoLib::AABB<GeoLib::Point> aabb(pnts_list.begin(), pnts_list.end());
 
@@ -87,14 +83,14 @@ TEST(GeoLibAABB, RandomNumberOfPointsRandomPointInAList)
 	 ASSERT_GE(half_box_size, max_pnt[2]) << "coordinate 2 of max_pnt is greater than " << half_box_size;
 }
 
-TEST(GeoLibAABB, RandomNumberOfPointsPointersToRandomPointsInAVector)
+TEST(GeoLibAABB, RandomNumberOfPointersToRandomPointsInAVector)
 {
 	/* initialize random seed: */
 	 srand ( time(NULL) );
 	 int n (rand() % 100000);
 	 int box_size (rand());
 	 int half_box_size(box_size/2);
-	 int minus_half_box_size(-1.0 * half_box_size);
+	 int minus_half_box_size(-half_box_size);
 
 	 // fill list with points
 	 std::vector<GeoLib::Point*> pnts;
@@ -102,8 +98,6 @@ TEST(GeoLibAABB, RandomNumberOfPointsPointersToRandomPointsInAVector)
 		 pnts.push_back(new GeoLib::Point(rand() % box_size - half_box_size, rand() % box_size - half_box_size, rand() % box_size - half_box_size));
 	 }
 
-	 std::cout << "testing with " << n << " points" << std::endl;
-
 	 // construct from list points a axis algined bounding box
 	 GeoLib::AABB<GeoLib::Point> aabb(pnts.begin(), pnts.end());
 
@@ -129,7 +123,7 @@ TEST(GeoLibAABB, RandomNumberOfPointsRandomPointInAVector)
 	 int n (rand() % 1000000);
 	 int box_size (rand());
 	 int half_box_size(box_size/2);
-	 int minus_half_box_size(-1.0 * half_box_size);
+	 int minus_half_box_size(-half_box_size);
 
 	 // fill list with points
 	 std::list<GeoLib::Point> pnts;
@@ -137,8 +131,6 @@ TEST(GeoLibAABB, RandomNumberOfPointsRandomPointInAVector)
 		 pnts.push_back(GeoLib::Point(rand() % box_size - half_box_size, rand() % box_size - half_box_size, rand() % box_size - half_box_size));
 	 }
 
-	 std::cout << "testing with " << n << " points" << std::endl;
-
 	 // construct from list points a axis algined bounding box
 	 GeoLib::AABB<GeoLib::Point> aabb(pnts.begin(), pnts.end());
 
@@ -153,3 +145,37 @@ TEST(GeoLibAABB, RandomNumberOfPointsRandomPointInAVector)
 	 ASSERT_GE(half_box_size, max_pnt[2]) << "coordinate 2 of max_pnt is greater than " << half_box_size;
 }
 
+TEST(GeoLibAABB, RandomNumberOfPointsRandomBox)
+{
+	/* initialize random seed: */
+	 srand (time(NULL));
+	 int n (rand() % 1000000);
+	 int box_size_x (rand());
+	 int box_size_y (rand());
+	 int box_size_z (rand());
+	 int half_box_size_x(box_size_x/2);
+	 int half_box_size_y(box_size_y/2);
+	 int half_box_size_z(box_size_z/2);
+	 int minus_half_box_size_x(-half_box_size_x);
+	 int minus_half_box_size_y(-half_box_size_y);
+	 int minus_half_box_size_z(-half_box_size_z);
+
+	 // fill list with points
+	 std::list<GeoLib::Point> pnts;
+	 for (int k(0); k<n; k++) {
+		 pnts.push_back(GeoLib::Point(rand() % box_size_x - half_box_size_x, rand() % box_size_y - half_box_size_y, rand() % box_size_z - half_box_size_z));
+	 }
+
+	 // construct from list points a axis algined bounding box
+	 GeoLib::AABB<GeoLib::Point> aabb(pnts.begin(), pnts.end());
+
+	 GeoLib::Point const& min_pnt(aabb.getMinPoint());
+	 GeoLib::Point const& max_pnt(aabb.getMaxPoint());
+
+	 ASSERT_LE(minus_half_box_size_x, min_pnt[0]) << "coordinate 0 of min_pnt is smaller than " << minus_half_box_size_x;
+	 ASSERT_LE(minus_half_box_size_y, min_pnt[1]) << "coordinate 1 of min_pnt is smaller than " << minus_half_box_size_y;
+	 ASSERT_LE(minus_half_box_size_z, min_pnt[2]) << "coordinate 2 of min_pnt is smaller than " << minus_half_box_size_z;
+	 ASSERT_GE(half_box_size_x, max_pnt[0]) << "coordinate 0 of max_pnt is greater than " << half_box_size_x;
+	 ASSERT_GE(half_box_size_y, max_pnt[1]) << "coordinate 1 of max_pnt is greater than " << half_box_size_y;
+	 ASSERT_GE(half_box_size_z, max_pnt[2]) << "coordinate 2 of max_pnt is greater than " << half_box_size_z;
+}