diff --git a/GeoLib/PointVec.cpp b/GeoLib/PointVec.cpp
index a9548aac5e44c7f35e1234c1b9704745f06bd236..c2cadc8812d5b4e67417204ebf33341238359f0f 100644
--- a/GeoLib/PointVec.cpp
+++ b/GeoLib/PointVec.cpp
@@ -74,30 +74,34 @@ void PointVec::push_back (Point* pnt, std::string const*const name)
 
 std::size_t PointVec::uniqueInsert (Point* pnt)
 {
-	std::size_t n (_data_vec->size()), k;
 	const double eps (std::numeric_limits<double>::epsilon());
-	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) {
-		_data_vec->push_back (pnt);
-		// update bounding box
-		_aabb.update (*((*_data_vec)[n]));
-		// update shortest distance
-		for (std::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;
-		}
-		return n;
+	auto const it = std::find_if(_data_vec->begin(), _data_vec->end(),
+		[&eps, &pnt](Point* const p)
+		{
+			return MathLib::maxNormDist(p, pnt) <= eps;
+		});
+
+	if (it != _data_vec->end())
+	{
+		delete pnt;
+		pnt = NULL;
+		return std::distance(_data_vec->begin(), it);
 	}
 
-	delete pnt;
-	pnt = NULL;
-	return k;
+	_data_vec->push_back(pnt);
+
+	// update bounding box
+	_aabb.update (*(_data_vec->back()));
+
+	// update shortest distance
+	std::for_each(_data_vec->begin(), _data_vec->end(),
+			[this](Point* const p)
+			{
+				_sqr_shortest_dist = std::min(_sqr_shortest_dist,
+						MathLib::sqrDist(p, _data_vec->back()));
+			});
+
+	return _data_vec->size()-1;
 }
 
 std::vector<Point*>* PointVec::filterStations(const std::vector<PropertyBounds> &bounds) const
@@ -134,10 +138,10 @@ void PointVec::makePntsUnique (std::vector<GeoLib::Point*>* pnt_vec,
 	// determine intervals with identical points to resort for stability of sorting
 	std::vector<std::size_t> identical_pnts_interval;
 	bool identical(false);
-	for (std::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 (std::size_t k = 0; k < n_pnts_in_file - 1; k++)
+	{
+		if (MathLib::maxNormDist((*pnt_vec)[k + 1], (*pnt_vec)[k]) <= eps)
+		{
 			// points are identical, sort by id
 			if (!identical)
 				identical_pnts_interval.push_back(k);
@@ -163,13 +167,9 @@ void PointVec::makePntsUnique (std::vector<GeoLib::Point*>* pnt_vec,
 	}
 
 	// check if there are identical points
-	for (std::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 (std::size_t k = 0; k < n_pnts_in_file - 1; k++)
+		if (MathLib::maxNormDist((*pnt_vec)[k + 1], (*pnt_vec)[k]) <= eps)
 			pnt_id_map[perm[k + 1]] = pnt_id_map[perm[k]];
-		}
-	}
 
 	// reverse permutation
 	BaseLib::Quicksort<std::size_t, GeoLib::Point*>(perm, 0, n_pnts_in_file, *pnt_vec);
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index abc76e3ebbf5fc730006c7d3c153c3d848b74734..15adac659646e539415fbf1e7fa65142db2984a7 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -118,6 +118,17 @@ T sqrDist(const MathLib::TemplatePoint<T>* p0, const MathLib::TemplatePoint<T>*
 /** squared dist between double arrays p0 and p1 (size of arrays is 3) */
 double sqrDist(const double* p0, const double* p1);
 
+/** Distance between points p0 and p1 in the maximum norm. */
+template <typename T>
+T maxNormDist(const MathLib::TemplatePoint<T>* p0, const MathLib::TemplatePoint<T>* p1)
+{
+	const T x = fabs((*p1)[0] - (*p0)[0]);
+	const T y = fabs((*p1)[1] - (*p0)[1]);
+	const T z = fabs((*p1)[2] - (*p0)[2]);
+
+	return std::max(x, std::max(y, z));
+}
+
 /** linear normalisation of val from [min, max] into [0,1] */
 float normalize(float min, float max, float val);
 
diff --git a/Tests/GeoLib/TestPointVec.cpp b/Tests/GeoLib/TestPointVec.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ae15d0e8993fd526ac9cfc75b2753cb49e85b7be
--- /dev/null
+++ b/Tests/GeoLib/TestPointVec.cpp
@@ -0,0 +1,132 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include "gtest/gtest.h"
+#include <ctime>
+#include <random>
+
+#include "GeoLib/PointVec.h"
+
+class PointVecTest : public testing::Test
+{
+public:
+	typedef std::vector<GeoLib::Point*> VectorOfPoints;
+
+	PointVecTest()
+		: gen(std::random_device() ()), name("JustAName")
+	{
+		ps_ptr = new VectorOfPoints;
+	}
+
+protected:
+	// Generates n new points according to given random number distribution,
+	// which is uniform distribution in [-1, 1]^3.
+	template <typename RandomDistribution = std::uniform_real_distribution<double>>
+	void
+	generateRandomPoints(std::size_t const n = 1000,
+		RandomDistribution rnd = std::uniform_real_distribution<double>(-1, 1))
+	{
+		std::generate_n(std::back_inserter(*ps_ptr), n,
+			[&]() { return new GeoLib::Point(rnd(gen), rnd(gen), rnd(gen)); });
+	}
+
+protected:
+	std::mt19937 gen;
+	const std::string name;
+
+	VectorOfPoints* ps_ptr;
+};
+
+// Testing nullptr input vector.
+TEST_F(PointVecTest, TestPointVecCtorNullptr)
+{
+	ASSERT_ANY_THROW(GeoLib::PointVec(name, nullptr));
+	delete ps_ptr;
+}
+
+// Testing empty input vector.
+TEST_F(PointVecTest, TestPointVecCtorEmpty)
+{
+	ASSERT_ANY_THROW(GeoLib::PointVec(name, ps_ptr));
+}
+
+// Testing input vector with single point.
+TEST_F(PointVecTest, TestPointVecCtorSinglePoint)
+{
+	ps_ptr->push_back(new GeoLib::Point(0,0,0));
+	GeoLib::PointVec* point_vec = nullptr;
+	ASSERT_NO_THROW(point_vec = new GeoLib::PointVec(name, ps_ptr));
+	ASSERT_EQ(std::size_t(1), point_vec->size());
+
+	delete point_vec;
+}
+
+// Testing input vector with two different points.
+TEST_F(PointVecTest, TestPointVecCtorTwoDiffPoints)
+{
+	ps_ptr->push_back(new GeoLib::Point(0,0,0));
+	ps_ptr->push_back(new GeoLib::Point(1,0,0));
+
+	GeoLib::PointVec* point_vec = nullptr;
+	ASSERT_NO_THROW(point_vec = new GeoLib::PointVec(name, ps_ptr));
+	ASSERT_EQ(std::size_t(2), point_vec->size());
+
+	delete point_vec;
+}
+
+// Testing input vector with two equal points.
+TEST_F(PointVecTest, TestPointVecCtorTwoEqualPoints)
+{
+	ps_ptr->push_back(new GeoLib::Point(0,0,0));
+	ps_ptr->push_back(new GeoLib::Point(0,0,0));
+
+	GeoLib::PointVec* point_vec = nullptr;
+	ASSERT_NO_THROW(point_vec = new GeoLib::PointVec(name, ps_ptr));
+	ASSERT_EQ(std::size_t(1), point_vec->size());
+
+	delete point_vec;
+}
+
+// Testing input vector with single point.
+TEST_F(PointVecTest, TestPointVecPushBack)
+{
+	ps_ptr->push_back(new GeoLib::Point(0,0,0));
+	ps_ptr->push_back(new GeoLib::Point(1,0,0));
+	ps_ptr->push_back(new GeoLib::Point(0,1,0));
+	ps_ptr->push_back(new GeoLib::Point(0,0,1));
+	GeoLib::PointVec point_vec(name, ps_ptr);
+
+	// Adding a new point with same coordinates changes nothing.
+	point_vec.push_back(new GeoLib::Point(0,0,0));
+	point_vec.push_back(new GeoLib::Point(1,0,0));
+	point_vec.push_back(new GeoLib::Point(0,1,0));
+	point_vec.push_back(new GeoLib::Point(0,0,1));
+
+	ASSERT_EQ(std::size_t(4), point_vec.size());
+}
+
+// Testing random input points.
+TEST_F(PointVecTest, TestPointVecCtorRandomPoints)
+{
+	generateRandomPoints(10000);
+
+	GeoLib::PointVec* point_vec = nullptr;
+	ASSERT_NO_THROW(point_vec = new GeoLib::PointVec(name, ps_ptr));
+
+	delete point_vec;
+}
+TEST_F(PointVecTest, TestPointVecCtorRandomPointsLargeEps)
+{
+	generateRandomPoints(10000);
+
+	GeoLib::PointVec* point_vec = nullptr;
+	ASSERT_NO_THROW(point_vec = new GeoLib::PointVec(name, ps_ptr,
+					nullptr, GeoLib::PointVec::POINT, 1e-2));
+
+	delete point_vec;
+}