diff --git a/GeoLib/AxisAlignedBoundingBox.cpp b/GeoLib/AxisAlignedBoundingBox.cpp
index b61a7c86d54f8611afffe3dc3c356150bf040537..7e42d1464b19ea1375f1cee57d4be1defca338f4 100644
--- a/GeoLib/AxisAlignedBoundingBox.cpp
+++ b/GeoLib/AxisAlignedBoundingBox.cpp
@@ -5,28 +5,31 @@
  *      Author: TF
  */
 
-#include <limits>
-#include <cstddef>
-#include <cmath>
 #include "AxisAlignedBoundingBox.h"
+#include <cmath>
+#include <cstddef>
+#include <limits>
 
-namespace GeoLib {
-
+namespace GeoLib
+{
 AABB::AABB ()
 {
-	for (std::size_t k(0); k<3; k++) {
+	for (std::size_t k(0); k < 3; k++)
+	{
 		_min_pnt[k] = std::numeric_limits<double>::max();
 		_max_pnt[k] = std::numeric_limits<double>::min();
 	}
 }
 
-AABB::AABB ( const std::vector<GeoLib::Point*> *points )
+AABB::AABB(AABB const& src) :
+	_min_pnt(src._min_pnt.getCoords()), _max_pnt(src._max_pnt.getCoords())
+{}
+
+AABB::AABB ( const std::vector<GeoLib::Point*>* points )
 {
 	size_t nPoints (points->size());
-	for (size_t i=0; i<nPoints; i++)
-	{
+	for (size_t i = 0; i < nPoints; i++)
 		this->update((*(*points)[i])[0], (*(*points)[i])[1], (*(*points)[i])[2]);
-	}
 }
 
 void AABB::update (GeoLib::Point const & pnt)
@@ -36,12 +39,18 @@ void AABB::update (GeoLib::Point const & pnt)
 
 void AABB::update (double x, double y, double z)
 {
-	if (x < _min_pnt[0]) _min_pnt[0] = x;
-	if (_max_pnt[0] < x) _max_pnt[0] = x;
-	if (y < _min_pnt[1]) _min_pnt[1] = y;
-	if (_max_pnt[1] < y) _max_pnt[1] = y;
-	if (z < _min_pnt[2]) _min_pnt[2] = z;
-	if (_max_pnt[2] < z) _max_pnt[2] = z;
+	if (x < _min_pnt[0])
+		_min_pnt[0] = x;
+	if (_max_pnt[0] < x)
+		_max_pnt[0] = x;
+	if (y < _min_pnt[1])
+		_min_pnt[1] = y;
+	if (_max_pnt[1] < y)
+		_max_pnt[1] = y;
+	if (z < _min_pnt[2])
+		_min_pnt[2] = z;
+	if (_max_pnt[2] < z)
+		_max_pnt[2] = z;
 }
 
 bool AABB::containsPoint (GeoLib::Point const & pnt, double eps) const
@@ -49,26 +58,27 @@ bool AABB::containsPoint (GeoLib::Point const & pnt, double eps) const
 	return containsPoint (pnt[0], pnt[1], pnt[2], eps);
 }
 
-bool AABB::containsPoint (const double *pnt, double eps) const
+bool AABB::containsPoint (const double* pnt, double eps) const
 {
 	return containsPoint (pnt[0], pnt[1], pnt[2], eps);
 }
 
 bool AABB::containsPoint (double x, double y, double z, double eps) const
 {
-	if ((_min_pnt[0] <= x && x <= _max_pnt[0])
-			|| std::fabs(_min_pnt[0]-x) < eps
-			|| std::fabs(x - _max_pnt[0]) < eps) {
-		if ((_min_pnt[1] <= y && y <= _max_pnt[1])
-				|| std::fabs(_min_pnt[1]-y) < eps
-				|| std::fabs(y-_max_pnt[1]) < eps) {
-			if ((_min_pnt[2] <= z && z <= _max_pnt[2])
-					|| std::fabs(_min_pnt[2]-z) < eps
-					|| std::fabs(z-_max_pnt[2]) < eps) {
+	if ((_min_pnt[0] <= x && x <= _max_pnt[0]) || std::fabs(_min_pnt[0] - x) < eps
+					|| std::fabs(x - _max_pnt[0]) < eps) {
+		if ((_min_pnt[1] <= y && y <= _max_pnt[1]) || std::fabs(_min_pnt[1] - y) < eps
+						|| std::fabs(y - _max_pnt[1]) < eps) {
+			if ((_min_pnt[2] <= z && z <= _max_pnt[2]) || std::fabs(_min_pnt[2] - z) < eps
+							|| std::fabs(z - _max_pnt[2]) < eps) {
 				return true;
-			} else return false;
-		} else return false;
+			} else {
+				return false;
+			}
+		} else {
+			return false;
+		}
 	} else return false;
 }
 
-}
+} // end namespace GeoLib
diff --git a/GeoLib/AxisAlignedBoundingBox.h b/GeoLib/AxisAlignedBoundingBox.h
index 1c0996fa4019f24d46b660853b35592b0a345917..844895fc0c20c1ba0481efcbcf1d9664d790dcec 100644
--- a/GeoLib/AxisAlignedBoundingBox.h
+++ b/GeoLib/AxisAlignedBoundingBox.h
@@ -9,11 +9,11 @@
 #define AXISALIGNEDBOUNDINGBOX_H_
 
 #include "Point.h"
-#include <vector>
 #include <limits>
+#include <vector>
 
-namespace GeoLib {
-
+namespace GeoLib
+{
 /**
  *
  * \ingroup GeoLib
@@ -28,10 +28,17 @@ public:
 	 * */
 	AABB ();
 
+	/**
+	 * copy constructor.
+	 * @param src an axis aligned bounding box
+	 * @return
+	 */
+	AABB(AABB const& src);
+
 	/**
 	 * construction of object using vector of points
 	 * */
-	AABB ( const std::vector<GeoLib::Point*> *points );
+	AABB ( const std::vector<GeoLib::Point*>* points );
 
 	void update (GeoLib::Point const & pnt);
 	/**
@@ -42,7 +49,7 @@ public:
 	/**
 	 * update axis aligned bounding box
 	 */
-	void update (const double *pnt)
+	void update (const double* pnt)
 	{
 		update (pnt[0], pnt[1], pnt[2]);
 	}
@@ -51,27 +58,29 @@ public:
 	 * check if point is in the axis aligned bounding box
 	 * (employing containsPoint (double x, double y, double z))
 	 */
-	bool containsPoint (GeoLib::Point const & pnt, double eps = std::numeric_limits<double>::epsilon()) const;
+	bool containsPoint (GeoLib::Point const & pnt,
+	                    double eps = std::numeric_limits<double>::epsilon()) const;
 
 	/**
 	 * wrapper for GeoLib::Point
 	 */
-	bool containsPoint (const double *pnt, double eps = std::numeric_limits<double>::epsilon()) const;
+	bool containsPoint (const double* pnt, double eps =
+	                            std::numeric_limits<double>::epsilon()) const;
 
 	/**
 	 * check if point described by its coordinates x, y, z is in
 	 * the axis aligned bounding box
 	 */
-	bool containsPoint (double x, double y, double z, double eps = std::numeric_limits<double>::epsilon()) const;
+	bool containsPoint(double x, double y, double z, double eps =
+					std::numeric_limits<double>::epsilon()) const;
 
-	GeoLib::Point getMinPoint () const { return _min_pnt; }
-	GeoLib::Point getMaxPoint () const { return _max_pnt; }
+	GeoLib::Point const& getMinPoint () const { return _min_pnt; }
+	GeoLib::Point const& getMaxPoint () const { return _max_pnt; }
 
-private:
+protected:
 	GeoLib::Point _min_pnt;
 	GeoLib::Point _max_pnt;
 };
-
 } // end namespace
 
 #endif /* AXISALIGNEDBOUNDINGBOX_H_ */
diff --git a/GeoLib/GEOObjects.cpp b/GeoLib/GEOObjects.cpp
index 866172c27aa978fd0e2f7b0c8ba574925d556bac..aa0aee8df3c36567a3967f558b4cc8908714b98d 100644
--- a/GeoLib/GEOObjects.cpp
+++ b/GeoLib/GEOObjects.cpp
@@ -5,7 +5,10 @@
  *      Author: TF / KR
  */
 
+// GeoLib
 #include "GEOObjects.h"
+
+// BaseLib
 #include "StringTools.h"
 
 #include <fstream>
@@ -500,7 +503,7 @@ int GEOObjects::exists(const std::string &geometry_name) const
 			return i;
 
 	// HACK for enabling conversion of files without loading the associated geometry
-	if (size>0 && _pnt_vecs[0]->getName().compare("conversionTestRun#1")==0)	
+	if (size>0 && _pnt_vecs[0]->getName().compare("conversionTestRun#1")==0)
 		return 1;
 
 	return -1;
diff --git a/GeoLib/PointVec.cpp b/GeoLib/PointVec.cpp
index e9a50bf93d1181b720000028ec30351b6c9ba811..143c107bb1e15a0144694835ce670adf11fb9373 100644
--- a/GeoLib/PointVec.cpp
+++ b/GeoLib/PointVec.cpp
@@ -5,46 +5,50 @@
  *      Author: TF
  */
 
-#include "PointVec.h"
+// GeoLib
 #include "BruteForceClosestPair.h"
+#include "PointVec.h"
+#include "PointWithID.h"
 
-namespace GeoLib {
+// MathLib
+#include "MathTools.h"
 
+namespace GeoLib
+{
 PointVec::PointVec (const std::string& name, std::vector<Point*>* points,
-		std::map<std::string, size_t>* name_id_map, PointType type) :
-			_pnt_vec (points), _name_id_map (name_id_map), _type(type), _name (name),
-			_sqr_shortest_dist (std::numeric_limits<double>::max())
+                    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())
 {
-	assert (_pnt_vec);
-	size_t number_of_all_input_pnts (_pnt_vec->size());
-
-	makePntsUnique (_pnt_vec, _pnt_id_map);
+	assert (_data_vec);
+	size_t number_of_all_input_pnts (_data_vec->size());
 
-	if (number_of_all_input_pnts - _pnt_vec->size() > 0)
-		std::cerr << "WARNING: there are " << number_of_all_input_pnts - _pnt_vec->size() << " double points" << std::endl;
-//	calculateShortestDistance ();
 	calculateAxisAlignedBoundingBox ();
+	rel_eps *= sqrt(MathLib::sqrDist (&(_aabb.getMinPoint()),&(_aabb.getMaxPoint())));
+	makePntsUnique (_data_vec, _pnt_id_map, rel_eps);
+
+	if (number_of_all_input_pnts - _data_vec->size() > 0)
+		std::cerr << "WARNING: there are " << number_of_all_input_pnts -
+		_data_vec->size() << " double points" << std::endl;
 }
 
 PointVec::~PointVec ()
-{
-	for (size_t k(0); k<_pnt_vec->size(); k++) {
-		delete (*_pnt_vec)[k];
-		(*_pnt_vec)[k] = NULL;
-	}
-	delete _pnt_vec;
-	delete _name_id_map;
-}
+{}
 
-size_t PointVec::push_back (Point *pnt)
+size_t PointVec::push_back (Point* pnt)
 {
 	_pnt_id_map.push_back (uniqueInsert(pnt));
-	return _pnt_id_map[_pnt_id_map.size()-1];
+	return _pnt_id_map[_pnt_id_map.size() - 1];
 }
 
-void PointVec::push_back (Point *pnt, const std::string& name)
+void PointVec::push_back (Point* pnt, std::string const*const name)
 {
-	std::map<std::string,size_t>::const_iterator it (_name_id_map->find (name));
+	if (name == NULL) {
+		_pnt_id_map.push_back (uniqueInsert(pnt));
+		return;
+	}
+
+	std::map<std::string,size_t>::const_iterator it (_name_id_map->find (*name));
 	if (it != _name_id_map->end()) {
 		std::cerr << "ERROR: PointVec::push_back (): two points with the same name" << std::endl;
 		return;
@@ -52,29 +56,28 @@ void PointVec::push_back (Point *pnt, const std::string& name)
 
 	size_t id (uniqueInsert (pnt));
 	_pnt_id_map.push_back (id);
-	(*_name_id_map)[name] = id;
+	(*_name_id_map)[*name] = id;
 }
 
 size_t PointVec::uniqueInsert (Point* pnt)
 {
-	size_t n (_pnt_vec->size()), k;
+	size_t n (_data_vec->size()), k;
 	const double eps (std::numeric_limits<double>::epsilon());
-	for (k=0; k<n; k++) {
-		if (fabs((*((*_pnt_vec)[k]))[0]-(*pnt)[0]) < eps
-				&&  fabs( (*((*_pnt_vec)[k]))[1]-(*pnt)[1]) < eps
-				&&  fabs( (*((*_pnt_vec)[k]))[2]-(*pnt)[2]) < eps) {
+	for (k = 0; k < n; k++)
+		if (fabs((*((*_data_vec)[k]))[0] - (*pnt)[0]) < eps
+		    &&  fabs( (*((*_data_vec)[k]))[1] - (*pnt)[1]) < eps
+		    &&  fabs( (*((*_data_vec)[k]))[2] - (*pnt)[2]) < eps)
 			break;
-		}
-	}
 
-	if(k==n) {
-		_pnt_vec->push_back (pnt);
-		// update shortest distance and bounding box
-		for (size_t i(0); i<n; i++) {
-			double sqr_dist (MathLib::sqrDist((*_pnt_vec)[i], (*_pnt_vec)[n]));
+	if(k == n) {
+		_data_vec->push_back (pnt);
+		// update bounding box
+		_aabb.update (*((*_data_vec)[n]));
+		// update shortest distance
+		for (size_t i(0); i < n; i++) {
+			double sqr_dist (MathLib::sqrDist((*_data_vec)[i], (*_data_vec)[n]));
 			if (sqr_dist < _sqr_shortest_dist)
 				_sqr_shortest_dist = sqr_dist;
-			aabb.update (*((*_pnt_vec)[n]));
 		}
 		return n;
 	}
@@ -84,73 +87,29 @@ size_t PointVec::uniqueInsert (Point* pnt)
 	return k;
 }
 
-std::vector<Point*> * PointVec::filterStations(const std::vector<PropertyBounds> &bounds) const
+std::vector<Point*>* PointVec::filterStations(const std::vector<PropertyBounds> &bounds) const
 {
-	std::vector<Point*> *tmpStations (new std::vector<Point*>);
-	size_t size (_pnt_vec->size());
-	for (size_t i=0; i<size; i++) {
-		if (static_cast<Station*>((*_pnt_vec)[i])->inSelection(bounds)) tmpStations->push_back((*_pnt_vec)[i]);
-	}
+	std::vector<Point*>* tmpStations (new std::vector<Point*>);
+	size_t size (_data_vec->size());
+	for (size_t i = 0; i < size; i++)
+		if (static_cast<Station*>((*_data_vec)[i])->inSelection(bounds))
+			tmpStations->push_back((*_data_vec)[i]);
 	return tmpStations;
 }
 
-bool PointVec::getElementIDByName (const std::string& name, size_t &id) const
-{
-	std::map<std::string,size_t>::const_iterator it (_name_id_map->find (name));
-
-	if (it != _name_id_map->end()) {
-		id = it->second;
-		return true;
-	} else return false;
-}
-
-const Point* PointVec::getElementByName (const std::string& name) const
-{
-	size_t id;
-	bool ret (getElementIDByName (name, id));
-	if (ret) {
-		return (*_pnt_vec)[id];
-	} else {
-		return NULL;
-	}
-}
-
-bool PointVec::getNameOfElement (const Point* data, std::string& name) const
-{
-	for (size_t k(0); k<_pnt_vec->size(); k++) {
-		if ((*_pnt_vec)[k] == data) {
-			return getNameOfElementByID (k, name);
-		}
-	}
-	return false;
-}
-
-bool PointVec::getNameOfElementByID (size_t id, std::string& element_name) const
-{
-	if (! _name_id_map) return false;
-	// search in map for id
-	std::map<std::string,size_t>::const_iterator it (_name_id_map->begin());
-	while (it != _name_id_map->end()) {
-		if (it->second == id) {
-			element_name = it->first;
-			return true;
-		}
-		it++;
-	}
-	return false;
-}
-
 double PointVec::getShortestPointDistance () const
 {
 	return sqrt (_sqr_shortest_dist);
 }
 
-void PointVec::makePntsUnique (std::vector<GeoLib::Point*>* pnt_vec, std::vector<size_t> &pnt_id_map)
+void PointVec::makePntsUnique (std::vector<GeoLib::Point*>* pnt_vec,
+                               std::vector<size_t> &pnt_id_map, double eps)
 {
 	size_t n_pnts_in_file (pnt_vec->size());
 	std::vector<size_t> perm;
 	pnt_id_map.reserve (n_pnts_in_file);
-	for (size_t k(0); k<n_pnts_in_file; k++) {
+	for (size_t k(0); k < n_pnts_in_file; k++)
+	{
 		perm.push_back (k);
 		pnt_id_map.push_back(k);
 	}
@@ -160,88 +119,130 @@ void PointVec::makePntsUnique (std::vector<GeoLib::Point*>* pnt_vec, std::vector
 
 	// unfortunately quicksort is not stable -
 	// sort identical points by id - to make sorting stable
-	double eps (sqrt(std::numeric_limits<double>::min()));
 	// determine intervals with identical points to resort for stability of sorting
 	std::vector<size_t> identical_pnts_interval;
 	bool identical (false);
-	for (size_t k=0; k<n_pnts_in_file-1; k++) {
-		if ( fabs((*((*pnt_vec)[k+1]))[0]-(*((*pnt_vec)[k]))[0]) < eps
-			&&  fabs( (*((*pnt_vec)[k+1]))[1]-(*((*pnt_vec)[k]))[1]) < eps
-			&&  fabs( (*((*pnt_vec)[k+1]))[2]-(*((*pnt_vec)[k]))[2]) < eps) {
+	for (size_t k = 0; k < n_pnts_in_file - 1; k++)
+	{
+		if ( fabs((*((*pnt_vec)[k + 1]))[0] - (*((*pnt_vec)[k]))[0]) < eps
+		     &&  fabs( (*((*pnt_vec)[k + 1]))[1] - (*((*pnt_vec)[k]))[1]) < eps
+		     &&  fabs( (*((*pnt_vec)[k + 1]))[2] - (*((*pnt_vec)[k]))[2]) < eps)
+		{
 			// points are identical, sort by id
-			if (!identical) identical_pnts_interval.push_back (k);
+			if (!identical)
+				identical_pnts_interval.push_back (k);
 			identical = true;
-		} else {
-			if (identical) identical_pnts_interval.push_back (k+1);
+		}
+		else
+		{
+			if (identical)
+				identical_pnts_interval.push_back (k + 1);
 			identical = false;
 		}
 	}
-	if (identical) identical_pnts_interval.push_back (n_pnts_in_file);
+	if (identical)
+		identical_pnts_interval.push_back (n_pnts_in_file);
 
-	for (size_t i(0); i<identical_pnts_interval.size()/2; i++) {
+	for (size_t i(0); i < identical_pnts_interval.size() / 2; i++)
+	{
 		// bubble sort by id
-		size_t beg (identical_pnts_interval[2*i]);
-		size_t end (identical_pnts_interval[2*i+1]);
-		for (size_t j (beg); j<end; j++) {
-			for (size_t k (beg); k<end-1; k++) {
-				if (perm[k] > perm[k+1]) std::swap (perm[k], perm[k+1]);
-			}
-		}
+		size_t beg (identical_pnts_interval[2 * i]);
+		size_t end (identical_pnts_interval[2 * i + 1]);
+		for (size_t j (beg); j < end; j++)
+			for (size_t k (beg); k < end - 1; k++)
+				if (perm[k] > perm[k + 1])
+					std::swap (perm[k], perm[k + 1]);
+
 	}
 
 	// check if there are identical points
-	for (size_t k=0; k<n_pnts_in_file-1; k++) {
-		if ( fabs((*((*pnt_vec)[k+1]))[0]-(*((*pnt_vec)[k]))[0]) < eps
-				&&  fabs( (*((*pnt_vec)[k+1]))[1]-(*((*pnt_vec)[k]))[1]) < eps
-				&&  fabs( (*((*pnt_vec)[k+1]))[2]-(*((*pnt_vec)[k]))[2]) < eps) {
-			pnt_id_map[perm[k+1]] = pnt_id_map[perm[k]];
-		}
-	}
+	for (size_t k = 0; k < n_pnts_in_file - 1; k++)
+		if ( fabs((*((*pnt_vec)[k + 1]))[0] - (*((*pnt_vec)[k]))[0]) < eps
+		     &&  fabs( (*((*pnt_vec)[k + 1]))[1] - (*((*pnt_vec)[k]))[1]) < eps
+		     &&  fabs( (*((*pnt_vec)[k + 1]))[2] - (*((*pnt_vec)[k]))[2]) < eps)
+			pnt_id_map[perm[k + 1]] = pnt_id_map[perm[k]];
 
 	// reverse permutation
 	BaseLib::Quicksort<GeoLib::Point*> (perm, 0, n_pnts_in_file, *pnt_vec);
 
 	// remove the second, third, ... occurrence from vector
-	for (size_t k(0); k<n_pnts_in_file; k++) {
-		if (pnt_id_map[k] < k) {
+	for (size_t k(0); k < n_pnts_in_file; k++)
+		if (pnt_id_map[k] < k)
+		{
 			delete (*pnt_vec)[k];
 			(*pnt_vec)[k] = NULL;
 		}
-	}
 	// remove NULL-ptr from vector
-	for (std::vector<GeoLib::Point*>::iterator it(pnt_vec->begin()); it != pnt_vec->end(); ) {
-		if (*it == NULL) {
+	for (std::vector<GeoLib::Point*>::iterator it(pnt_vec->begin()); it != pnt_vec->end(); )
+	{
+		if (*it == NULL)
 			it = pnt_vec->erase (it);
-		}
-		else it++;
+		else
+			it++;
 	}
 
 	// renumber id-mapping
 	size_t cnt (0);
-	for (size_t k(0); k<n_pnts_in_file; k++) {
-		if (pnt_id_map[k] == k) { // point not removed, if necessary: id change
+	for (size_t k(0); k < n_pnts_in_file; k++)
+	{
+		if (pnt_id_map[k] == k) // point not removed, if necessary: id change
+		{
 			pnt_id_map[k] = cnt;
 			cnt++;
-		} else {
-			pnt_id_map[k] = pnt_id_map[pnt_id_map[k]];
 		}
+		else
+			pnt_id_map[k] = pnt_id_map[pnt_id_map[k]];
 	}
 
+	// KR correct renumbering of indices
+//	size_t cnt(0);
+//	std::map<size_t, size_t> reg_ids;
+//	for (size_t k(0); k < n_pnts_in_file; k++) {
+//		if (pnt_id_map[k] == k) {
+//			reg_ids.insert(std::pair<size_t, size_t>(k, cnt));
+//			cnt++;
+//		} else reg_ids.insert(std::pair<size_t, size_t>(k, reg_ids[pnt_id_map[k]]));
+//	}
+//	for (size_t k(0); k < n_pnts_in_file; k++)
+//		pnt_id_map[k] = reg_ids[k];
 }
 
 void PointVec::calculateShortestDistance ()
 {
 	size_t i, j;
-	BruteForceClosestPair (*_pnt_vec, i, j);
-	_sqr_shortest_dist = MathLib::sqrDist ((*_pnt_vec)[i], (*_pnt_vec)[j]);
+	BruteForceClosestPair (*_data_vec, i, j);
+	_sqr_shortest_dist = MathLib::sqrDist ((*_data_vec)[i], (*_data_vec)[j]);
 }
 
 void PointVec::calculateAxisAlignedBoundingBox ()
 {
-	const size_t n_pnts (_pnt_vec->size());
-	for (size_t i(0); i<n_pnts; i++) {
-		aabb.update (*(*_pnt_vec)[i]);
-	}
+	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()));
+
+	const size_t nPoints(subset.size());
+	for (size_t i = 0; i < nPoints; i++)
+		(*new_points)[i] = new GeoLib::PointWithID((*this->_data_vec)[subset[i]]->getCoords(), subset[i]);
+
+	return new_points;
+}
+
+std::vector<GeoLib::Point*>* PointVec::deepcopy(const std::vector<GeoLib::Point*> *pnt_vec)
+{
+	std::vector<GeoLib::Point*>* new_points (new std::vector<GeoLib::Point*>);
+
+	const size_t nPoints(pnt_vec->size());
+	for (size_t i = 0; i < nPoints; i++)
+		new_points->push_back(new GeoLib::Point((*pnt_vec)[i]->getCoords()));
+	return new_points;
+}
+
+
+
+
 } // end namespace
diff --git a/GeoLib/PointVec.h b/GeoLib/PointVec.h
index 6e7282882b426314e6ef6a96dfd5b4cb53f3d1df..4d6c47f51e9056ed406c561264b667597003ff37 100644
--- a/GeoLib/PointVec.h
+++ b/GeoLib/PointVec.h
@@ -5,24 +5,28 @@
  *      Author: TF / KR
  */
 
-
 // GeoLib
+#include "AxisAlignedBoundingBox.h"
 #include "Point.h"
 #include "Station.h"
-#include "AxisAlignedBoundingBox.h"
 
-// Base
-#include "quicksort.h"
+// BaseLib
 #include "binarySearch.h"
+#include "quicksort.h"
 
-#include <vector>
-#include <string>
 #include <map>
+#include <string>
+#include <vector>
 
 #ifndef POINTVEC_H_
 #define POINTVEC_H_
 
-namespace GeoLib {
+#include "TemplateVec.h"
+
+namespace GeoLib
+{
+
+class PointWithID;
 
 /**
  * \ingroup GeoLib
@@ -32,7 +36,7 @@ namespace GeoLib {
  * a unique name from class GEOObject. For this reason PointVec should have
  * a name.
  * */
-class PointVec
+class PointVec : public TemplateVec<Point>
 {
 public:
 	/// Signals if the vector contains object of type Point or Station
@@ -55,13 +59,18 @@ public:
 	 * PointVec will take the ownership of the vector, i.e. it
 	 * deletes the names
 	 * @param type the type of the point, \sa enum PointType
+	 * @param rel_eps This is a relative error tolerance value for the test of identical points.
+	 * The size of the axis aligned bounding box multiplied with the value of rel_eps gives the
+	 * real tolerance \f$tol\f$. Two points \f$p_0, p_1 \f$ are identical iff
+	 * \f$|p_1 - p_0| \le tol.\f$
 	 * @return an object of type PointVec
 	 */
-	PointVec (const std::string& name, std::vector<Point*>* points, std::map<std::string, size_t>* name_id_map = NULL,
-				PointType type = PointVec::POINT);
+	PointVec (const std::string& name, std::vector<Point*>* points,
+	          std::map<std::string, size_t>* name_id_map = NULL,
+	          PointType type = PointVec::POINT, double rel_eps = sqrt(std::numeric_limits<double>::min()));
 
 	/** Destructor deletes all Points of this PointVec. */
-	~PointVec ();
+	virtual ~PointVec ();
 
 	/**
 	 * Method adds a Point to the (internal) standard vector and takes the ownership.
@@ -70,83 +79,36 @@ public:
 	 * @param pnt the pointer to the Point
 	 * @return the id of the point within the internal vector
 	 */
-	size_t push_back (Point *pnt);
+	size_t push_back (Point* pnt);
 
 	/**
-	 * push_back adds new elements at the end of the vector _pnt_vec.
+	 * push_back adds new elements at the end of the vector _data_vec.
 	 * @param pnt a pointer to the point, PointVec takes ownership of the point
 	 * @param name the name of the point
 	 */
-	void push_back (Point *pnt, const std::string& name);
+	virtual void push_back (Point* pnt, std::string const*const name);
 
-	/**
-	 * get the actual number of Points
-	 */
-	size_t size () const { return _pnt_vec->size(); }
 	/**
 	 * get the type of Point, this can be either POINT or STATION
 	 *
 	 */
 	PointType getType() const { return _type; }
 
-	/**
-	 * getVector returns the internal vector of Points,
-	 * you are not able to change the Points or the address of the vector.
-	 */
-	const std::vector<Point*>* getVector () const { return _pnt_vec; }
-
-	std::vector<Point*> *filterStations(const std::vector<PropertyBounds> &bounds) const;
-
-	/** sets the name of the object
-	 * \param n the name as standard string */
-	void setName(const std::string & n) { _name = n; }
-	/** returns the name of the object */
-	std::string getName () const { return _name; }
-
-	/**
-	 * search the vector of names for the ID of the point with the given name
-	 * @param name the name of the point
-	 * @param id the id of the point
-	 * @return the id of the point
-	 */
-	bool getElementIDByName (const std::string& name, size_t &id) const;
-
-	/**
-	 * Method searchs for point with the given name. If it found a point with the name
-	 * it returns a pointer to the point, else it returns the NULL pointer.
-	 * @param name the name of the point
-	 * @return the pointer to the point or NULL
-	 */
-	const Point* getElementByName (const std::string& name) const;
-
-	/**
-	 * The method returns true if the given element of type T
-	 * can be found and the element has a name, else method returns false.
-	 * @param data the data element, one wants to know the name
-	 * @param name the name of the data element (if the data element is
-	 * found and the data element has a name)
-	 * @return if element is found and has a name: true, else false
-	 */
-	bool getNameOfElement (const Point* data, std::string& name) const;
-
-	/**
-	 * The method returns true if there is a name associated
-	 * with the given id, else method returns false.
-	 * @param id the id
-	 * @param element_name if a name associated with the id
-	 * is found name is assigned to element_name
-	 * @return if there is name associated with the given id:
-	 * true, else false
-	 */
-	bool getNameOfElementByID (size_t id, std::string& element_name) const;
+	std::vector<Point*>* filterStations(const std::vector<PropertyBounds> &bounds) const;
 
 	const std::vector<size_t>& getIDMap () const { return _pnt_id_map; }
 
 	double getShortestPointDistance () const;
 	const GeoLib::AABB& getAxisAlignedBoundingBox () const;
 
+	/// Creates a real copy of the point vector in memeory.
+	static std::vector<GeoLib::Point*>* deepcopy(const std::vector<GeoLib::Point*> *pnt_vec);
+
+	/// Returns a subset of this point vector containing only the points specified in "subset" as PointWithID-objects
+	std::vector<GeoLib::Point*>* getSubset(const std::vector<size_t> &subset);
+
 private:
-	void makePntsUnique (std::vector<GeoLib::Point*>* pnt_vec, std::vector<size_t> &pnt_id_map);
+	void makePntsUnique (std::vector<GeoLib::Point*>* pnt_vec, std::vector<size_t> &pnt_id_map, double eps = sqrt(std::numeric_limits<double>::min()));
 
 	/** copy constructor doesn't have an implementation */
 	// compiler does not create a (possible unwanted) copy constructor
@@ -159,23 +121,11 @@ private:
 	// this way the compiler does not create a (possible unwanted) assignment operator
 	PointVec& operator= (const PointVec& rhs);
 
-	size_t uniqueInsert (Point *pnt);
+	size_t uniqueInsert (Point* pnt);
 
-	/**
-	 * pointer to a vector of pointers to Points
-	 *
-	 * The destructor of PointVec will delete all GeoLib::Points
-	 * inside the vector.
-	 */
-	std::vector<Point*> *_pnt_vec;
-	/**
-	 * used to store the name associated with a point
-	 */
-	std::map<std::string, size_t>* _name_id_map;
 	/** the type of the point (\sa enum PointType) */
 	PointType _type;
-	/** the name of the object */
-	std::string _name;
+
 	/**
 	 * permutation of the geometric elements according
 	 * to their lexicographical order
@@ -192,9 +142,8 @@ private:
 	double _sqr_shortest_dist;
 
 	void calculateAxisAlignedBoundingBox ();
-	AABB aabb;
+	AABB _aabb;
 };
-
 } // end namespace
 
 #endif /* POINTVEC_H_ */
diff --git a/GeoLib/TemplateVec.h b/GeoLib/TemplateVec.h
index 51ab39ef0b31562d584cb5588d22bd8821206d96..1406f6e5de7ca4bef165917e8aeffd84d88a5f33 100644
--- a/GeoLib/TemplateVec.h
+++ b/GeoLib/TemplateVec.h
@@ -8,8 +8,8 @@
 #ifndef TEMPLATEVEC_H_
 #define TEMPLATEVEC_H_
 
-namespace GeoLib {
-
+namespace GeoLib
+{
 /**
  * \ingroup GeoLib
  *
@@ -32,7 +32,8 @@ public:
 	 * the data_vec.
 
 	 */
-	TemplateVec (const std::string &name, std::vector<T*> *data_vec, std::map<std::string, size_t>* elem_name_map = NULL) :
+	TemplateVec (const std::string &name, std::vector<T*>* data_vec,
+	             std::map<std::string, size_t>* elem_name_map = NULL) :
 		_name(name), _data_vec(data_vec), _name_id_map (elem_name_map)
 	{}
 
@@ -41,7 +42,7 @@ public:
 	 */
 	virtual ~TemplateVec ()
 	{
-		for (size_t k(0); k<size(); k++) delete (*_data_vec)[k];
+		for (size_t k(0); k < size(); k++) delete (*_data_vec)[k];
 		delete _data_vec;
 		delete _name_id_map;
 	}
@@ -77,21 +78,23 @@ public:
 	{
 		std::map<std::string,size_t>::const_iterator it (_name_id_map->find (name));
 
-		if (it != _name_id_map->end()) {
+		if (it != _name_id_map->end())
+		{
 			id = it->second;
 			return true;
-		} else return false;
+		}
+		else return false;
 	}
 
+	/// Returns an element with the given name.
 	const T* getElementByName (const std::string& name) const
 	{
 		size_t id;
 		bool ret (getElementIDByName (name, id));
-		if (ret) {
+		if (ret)
 			return (*_data_vec)[id];
-		} else {
+		else
 			return NULL;
-		}
 	}
 
 	/**
@@ -105,11 +108,13 @@ public:
 	 */
 	bool getNameOfElementByID (size_t id, std::string& element_name) const
 	{
-		if (! _name_id_map) return false;
+		if (!_name_id_map) return false;
 		// search in map for id
 		std::map<std::string,size_t>::const_iterator it (_name_id_map->begin());
-		while (it != _name_id_map->end()) {
-			if (it->second == id) {
+		while (it != _name_id_map->end())
+		{
+			if (it->second == id)
+			{
 				element_name = it->first;
 				return true;
 			}
@@ -118,9 +123,10 @@ public:
 		return false;
 	}
 
+	/// Return the name of an element based on its ID.
 	void setNameOfElementByID (size_t id, std::string& element_name)
 	{
-		if (! _name_id_map) return;
+		if (!_name_id_map) return;
 		_name_id_map->insert( std::pair<std::string, size_t>(element_name, id) );
 	}
 
@@ -134,28 +140,48 @@ public:
 	 */
 	bool getNameOfElement (const T* data, std::string& name) const
 	{
-		for (size_t k(0); k<_data_vec->size(); k++) {
-			if ((*_data_vec)[k] == data) {
+		for (size_t k(0); k < _data_vec->size(); k++)
+			if ((*_data_vec)[k] == data)
 				return getNameOfElementByID (k, name);
-			}
-		}
+
 		return false;
 	}
 
-	void push_back (T* data_element, std::string const * const name = NULL)
+	/// Adds a new element to the vector.
+	virtual void push_back (T* data_element, std::string const* const name = NULL)
 	{
 		_data_vec->push_back (data_element);
 		if (name == NULL) return;
-		if (! name->empty()) {
-			if (_name_id_map == NULL) {
+		if (!name->empty())
+		{
+			if (_name_id_map == NULL)
 				_name_id_map = new std::map <std::string, size_t>;
-			}
-			_name_id_map->insert (std::pair<std::string,size_t>(*name, _data_vec->size()-1));
+			_name_id_map->insert (std::pair<std::string,size_t>(*name, _data_vec->size() - 1));
 		}
 	}
 
+	/// Sets the given name for the element of the given ID.
+	void setNameForElement(size_t id, std::string const& name)
+	{
+		if (_name_id_map == NULL)
+			_name_id_map = new std::map<std::string, size_t>;
+
+		if ( !_name_id_map->empty())
+		{
+			for (std::map<std::string, size_t>::iterator it = _name_id_map->begin(); it != _name_id_map->end(); ++it)
+				if (it->second == id)
+				{
+					_name_id_map->erase(it); //check if old name already exists and delete it
+					break;
+				}
+		}
+		if (!name.empty()) {
+			//insert new or revised name
+			_name_id_map->insert(std::pair<std::string, size_t>(name, id));
+		}
+	}
 
-private:
+protected:
 	/** copy constructor doesn't have an implementation */
 	// compiler does not create a (possible unwanted) copy constructor
 	TemplateVec (const TemplateVec &);
@@ -169,13 +195,12 @@ private:
 	/**
 	 * pointer to a vector of data elements
 	 */
-	std::vector <T*> *_data_vec;
+	std::vector <T*>* _data_vec;
 	/**
 	 * store names associated with the element ids
 	 */
 	std::map<std::string, size_t>* _name_id_map;
 };
-
 } // end namespace GeoLib
 
 #endif /* TEMPLATEVEC_H_ */