diff --git a/GeoLib/GEOObjects.cpp b/GeoLib/GEOObjects.cpp
index 6a25fd445cbbbf368a7b92b8c6d93a70b07fc3fd..afa1ac054d0738e8b9074fefd39e1bc1e89035bf 100644
--- a/GeoLib/GEOObjects.cpp
+++ b/GeoLib/GEOObjects.cpp
@@ -404,27 +404,32 @@ bool GEOObjects::mergePoints(std::vector<std::string> const & geo_names,
 	std::vector<GeoLib::Point*>* merged_points (new std::vector<GeoLib::Point*>);
 	std::map<std::string, std::size_t>* merged_pnt_names(new std::map<std::string, std::size_t>);
 
-	for (std::size_t j(0); j < n_geo_names; j++) {
-		const std::vector<GeoLib::Point*>* pnts(this->getPointVec(geo_names[j]));
-		if (pnts) {
-			std::size_t n_pnts(0);
-			// do not consider stations
-			if (!dynamic_cast<GeoLib::Station*>((*pnts)[0])) {
-				std::string tmp_name;
-				n_pnts = pnts->size();
-				for (std::size_t k(0); k < n_pnts; k++) {
-					merged_points->push_back(new GeoLib::Point(((*pnts)[k])->getCoords()));
-					if (this->getPointVecObj(geo_names[j])->getNameOfElementByID(k, tmp_name)) {
-						merged_pnt_names->insert(
-								std::pair<std::string, std::size_t>(tmp_name, pnt_offsets[j] + k));
-					}
-				}
-			}
-			if (n_geo_names - 1 > j) {
-				pnt_offsets[j + 1] = n_pnts + pnt_offsets[j];
+	for (std::size_t j(0); j < n_geo_names; ++j) {
+		GeoLib::PointVec const*const pnt_vec(this->getPointVecObj(geo_names[j]));
+		if (pnt_vec == nullptr)
+			continue;
+		const std::vector<GeoLib::Point*>* pnts(pnt_vec->getVector());
+		if (pnts == nullptr) {
+			return false;
+		}
+
+		// do not consider stations
+		if (dynamic_cast<GeoLib::Station*>((*pnts)[0])) {
+			continue;
+		}
+
+		std::size_t const n_pnts(pnts->size());
+		for (std::size_t k(0); k < n_pnts; ++k) {
+			merged_points->push_back(new GeoLib::Point(*(*pnts)[k]));
+			std::string const& item_name(pnt_vec->getItemNameByID(k));
+			if (! item_name.empty()) {
+				merged_pnt_names->insert(
+					std::make_pair(item_name, pnt_offsets[j] + k));
 			}
-		} else
-			return false; //if no points for a given geometry are found, something is fundamentally wrong
+		}
+		if (n_geo_names - 1 > j) {
+			pnt_offsets[j + 1] = n_pnts + pnt_offsets[j];
+		}
 	}
 
 	addPointVec (merged_points, merged_geo_name, merged_pnt_names, 1e-6);
diff --git a/GeoLib/PointVec.cpp b/GeoLib/PointVec.cpp
index d53951f00e7707d538d3a2ef1f8f595a657c406d..a65c76af4101b5b7732157bd43cf2f9ed1034117 100644
--- a/GeoLib/PointVec.cpp
+++ b/GeoLib/PointVec.cpp
@@ -32,13 +32,14 @@ PointVec::PointVec (const std::string& name, std::vector<Point*>* points,
                     std::map<std::string, std::size_t>* name_id_map, PointType type, double rel_eps) :
 	TemplateVec<Point> (name, points, name_id_map),
 	_type(type),
-	_aabb(points->begin(), points->end())
+	_aabb(points->begin(), points->end()),
+	_rel_eps(rel_eps)
 {
 	assert (_data_vec);
 	std::size_t const number_of_all_input_pnts (_data_vec->size());
 
-	rel_eps *= sqrt(MathLib::sqrDist (_aabb.getMinPoint(),_aabb.getMaxPoint()));
-	makePntsUnique (_data_vec, _pnt_id_map, rel_eps);
+	_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())
 		WARN("PointVec::PointVec(): there are %d double points.",
@@ -86,11 +87,10 @@ void PointVec::push_back (Point* pnt, std::string const*const name)
 
 std::size_t PointVec::uniqueInsert (Point* pnt)
 {
-	const double eps (std::numeric_limits<double>::epsilon());
 	auto const it = std::find_if(_data_vec->begin(), _data_vec->end(),
-		[&eps, &pnt](Point* const p)
+		[this, &pnt](Point* const p)
 		{
-			return MathLib::maxNormDist(p, pnt) <= eps;
+			return MathLib::maxNormDist(p, pnt) <= _rel_eps;
 		});
 
 	if (it != _data_vec->end())
diff --git a/GeoLib/PointVec.h b/GeoLib/PointVec.h
index 291ccf109fdd62981765fb3656dd3dc6f1293042..758ba49694eac8c7b0bf288c40ee80d5c4fa361d 100644
--- a/GeoLib/PointVec.h
+++ b/GeoLib/PointVec.h
@@ -145,6 +145,7 @@ private:
 	std::vector<std::string> _id_to_name_map;
 
 	AABB<GeoLib::Point> _aabb;
+	double _rel_eps;
 };
 } // end namespace