diff --git a/GeoLib/Point.cpp b/GeoLib/Point.cpp
index abdd9590de73bb0a0761eaaf4c708d2482a700ba..92efe2c8c87ff51313e03dc3599c480e6783039f 100644
--- a/GeoLib/Point.cpp
+++ b/GeoLib/Point.cpp
@@ -12,48 +12,66 @@
 
 
 #include <cmath>
-#include <limits>
 
 #include "Point.h"
 
 namespace GeoLib {
 
-bool lessX (GeoLib::Point const & p0, GeoLib::Point const & p1)
+bool operator<= (const GeoLib::Point& p0, const GeoLib::Point& p1)
 {
-	if (p0[0] <= p1[0]) return true;
-	return false;
-}
+	if (p0[0] > p1[0]) {
+		return false;
+	} else {
+		if (p0[0] < p1[0]) {
+			return true;
+		}
+	}
+	// => p0[0] == p1[0]
 
-bool lessY (GeoLib::Point const & p0, GeoLib::Point const & p1)
-{
-	if (p0[1] <= p1[1]) return true;
-	return false;
-}
+	if (p0[1] > p1[1]) {
+		return false;
+	} else {
+		if (p0[1] < p1[1]) {
+			return true;
+		}
+	}
+	// => p0[1] == p1[1]
 
-bool lessZ (GeoLib::Point const & p0, GeoLib::Point const & p1)
-{
-	if (p0[2] <= p1[2]) return true;
-	return false;
+	if (p0[2] > p1[2]) {
+		return false;
+	} else {
+		return true;
+	}
 }
 
-bool operator<= (const GeoLib::Point& p0, const GeoLib::Point& p1)
+bool lessEq(const GeoLib::Point& p0, const GeoLib::Point& p1, double tol)
 {
-	double tol (sqrt (std::numeric_limits<double>::min()));
-
-	if (fabs (p0[0]-p1[0]) > tol * fabs(p0[0])) {
-		if (p0[0] < p1[0]) return true;
-		else return false;
+	// test a relative and an absolute criterion
+	if (fabs(p0[0]-p1[0]) > tol * std::min(fabs(p1[0]), fabs(p0[0])) && fabs(p0[0]-p1[0]) > tol) {
+		if (p0[0] <= p1[0])
+			return true;
+		else
+			return false;
 	} else {
 		// assume p0[0] == p1[0]
-		if (fabs (p0[1]-p1[1]) > tol * fabs(p0[1])) {
-			if (p0[1] < p1[1]) return true;
-			else return false;
+		if (fabs (p0[1]-p1[1]) > tol * fabs(p0[1]) && fabs(p0[1]-p1[1]) > tol) {
+			if (p0[1] <= p1[1])
+				return true;
+			else
+				return false;
 		} else {
 			// assume p0[1] == p1[1] and p0[0] == p1[0]
-			if (p0[2] < p1[2]) return true;
-			else return false;
+			if (fabs (p0[2]-p1[2]) > tol * fabs(p0[2]) && fabs(p0[2]-p1[2]) > tol) {
+				if (p0[2] <= p1[2])
+					return true;
+				else
+					return false;
+			} else {
+				return true;
+			}
 		}
 	}
 }
 
+
 } // end namespace GeoLib
diff --git a/GeoLib/Point.h b/GeoLib/Point.h
index e793ef5b3e7fc818cf679838c210f9ee628d96a4..cd4c0728344085eff867eedf73fc5f9dc57bcf0a 100644
--- a/GeoLib/Point.h
+++ b/GeoLib/Point.h
@@ -13,6 +13,8 @@
 #ifndef POINT_H_
 #define POINT_H_
 
+#include <limits>
+
 #include "TemplatePoint.h"
 
 namespace GeoLib {
@@ -24,33 +26,19 @@ namespace GeoLib {
 typedef TemplatePoint<double> Point;
 
 /**
- * comparison based on the x coordinate
- * @param p0 first point
- * @param p1 second point
- * @return true if the x coordinate of p0 is smaller equal the x coordinate of p1, else false
- */
-bool lessX (Point const & p0, Point const & p1);
-
-/**
- * comparison based on the y coordinate
- * @param p0 first point
- * @param p1 second point
- * @return true if the y coordinate of p0 is smaller equal the y coordinate of p1, else false
- */
-bool lessY (Point const & p0, Point const & p1);
-
-/**
- * comparison based on the z coordinate
- * @param p0 first point
- * @param p1 second point
- * @return true if the z coordinate of p0 is smaller equal the z coordinate of p1, else false
+ * lexicographic comparison of points
  */
-bool lessZ (Point const & p0, Point const & p1);
+bool operator<= (GeoLib::Point const & p0, GeoLib::Point const & p1);
 
 /**
- * lexicographic comparison of points
+ * lexicographical comparison of points taking an epsilon into account
+ * @param p0 first input Point
+ * @param p1 first input Point
+ * @param tol tolerance (if in the comparison operation the property fabs(p0[k] - p1[k]) < tol
+ * 	holds for the k-th coordinate the points are assumed the be equal in this coordinate)
+ * @return true, if p0 is lexicographically smaller than p1
  */
-bool operator<= (GeoLib::Point const & p0, GeoLib::Point const & p1);
+bool lessEq(const GeoLib::Point& p0, const GeoLib::Point& p1, double tol = std::numeric_limits<double>::epsilon());
 }
 
 
diff --git a/Tests/GeoLib/TestPoint.cpp b/Tests/GeoLib/TestPoint.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ef0c08f3fca15462151f483f716adc7a687fa73a
--- /dev/null
+++ b/Tests/GeoLib/TestPoint.cpp
@@ -0,0 +1,157 @@
+/**
+ * @file TestPoint.cpp
+ * @author Thomas Fischer
+ * @date Nov 8, 2012
+ *
+ * @copyright
+ * Copyright (c) 2012, 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 <ctime>
+#include "gtest/gtest.h"
+
+#include "Point.h"
+
+using namespace GeoLib;
+
+TEST(GeoLib, PointComparisonLessEq)
+{
+	// first coordinate
+	ASSERT_TRUE(lessEq(Point(0,1,1),Point(1,1,1)));
+	ASSERT_FALSE(lessEq(Point(1,1,1),Point(0,1,1)));
+	// second coordinate
+	ASSERT_TRUE(lessEq(Point(1,0,1),Point(1,1,1)));
+	ASSERT_FALSE(lessEq(Point(1,1,1),Point(1,0,1)));
+	// third coordinate
+	ASSERT_TRUE(lessEq(Point(1,1,0),Point(1,1,1)));
+	ASSERT_FALSE(lessEq(Point(1,1,1),Point(1,1,0)));
+
+	const double e(2*std::numeric_limits<double>::epsilon());
+	// first coordinate
+	ASSERT_TRUE(lessEq(Point(1-e,1,1),Point(1,1,1)));
+	ASSERT_FALSE(lessEq(Point(1,1,1),Point(1-e,1,1)));
+	// second coordinate
+	ASSERT_TRUE(lessEq(Point(1,1-e,1),Point(1,1,1)));
+	ASSERT_FALSE(lessEq(Point(1,1,1),Point(1,1-e,1)));
+	// third coordinate
+	ASSERT_TRUE(lessEq(Point(1,1,1-e),Point(1,1,1)));
+	ASSERT_FALSE(lessEq(Point(1,1,1),Point(1,1,1-e)));
+
+	ASSERT_TRUE(lessEq(Point(1,1,1),Point(1,1,1)));
+
+	const double half_eps(0.5*std::numeric_limits<double>::epsilon());
+	ASSERT_TRUE(lessEq(Point(1+half_eps,1,1),Point(1,1,1)));
+	ASSERT_TRUE(lessEq(Point(1,1+half_eps,1),Point(1,1,1)));
+	ASSERT_TRUE(lessEq(Point(1,1,1+half_eps),Point(1,1,1)));
+	ASSERT_TRUE(lessEq(Point(1,1,1),Point(1+half_eps,1,1)));
+	ASSERT_TRUE(lessEq(Point(1,1,1),Point(1,1+half_eps,1)));
+	ASSERT_TRUE(lessEq(Point(1,1,1),Point(1,1,1+half_eps)));
+
+	ASSERT_TRUE(lessEq(Point(1-half_eps,1,1),Point(1,1,1)));
+	ASSERT_TRUE(lessEq(Point(1,1-half_eps,1),Point(1,1,1)));
+	ASSERT_TRUE(lessEq(Point(1,1,1-half_eps),Point(1,1,1)));
+	ASSERT_TRUE(lessEq(Point(1,1,1),Point(1-half_eps,1,1)));
+	ASSERT_TRUE(lessEq(Point(1,1,1),Point(1,1-half_eps,1)));
+	ASSERT_TRUE(lessEq(Point(1,1,1),Point(1,1,1-half_eps)));
+
+	const double m(std::numeric_limits<double>::min());
+	ASSERT_TRUE(lessEq(Point(m+half_eps,m,m),Point(m,m,m)));
+	ASSERT_TRUE(lessEq(Point(m,m+half_eps,m),Point(m,m,m)));
+	ASSERT_TRUE(lessEq(Point(m,m,m+half_eps),Point(m,m,m)));
+	ASSERT_TRUE(lessEq(Point(m,m,m),Point(m+half_eps,m,m)));
+	ASSERT_TRUE(lessEq(Point(m,m,m),Point(m,m+half_eps,m)));
+	ASSERT_TRUE(lessEq(Point(m,m,m),Point(m,m,m+half_eps)));
+
+	const double zero(0.0);
+	ASSERT_TRUE(lessEq(Point(zero+half_eps,zero,zero),Point(zero,zero,zero)));
+	ASSERT_TRUE(lessEq(Point(zero,zero+half_eps,zero),Point(zero,zero,zero)));
+	ASSERT_TRUE(lessEq(Point(zero,zero,zero+half_eps),Point(zero,zero,zero)));
+	ASSERT_TRUE(lessEq(Point(zero,zero,zero),Point(zero+half_eps,zero,zero)));
+	ASSERT_TRUE(lessEq(Point(zero,zero,zero),Point(zero,zero+half_eps,zero)));
+	ASSERT_TRUE(lessEq(Point(zero,zero,zero),Point(zero,zero,zero+half_eps)));
+
+	ASSERT_TRUE(lessEq(Point(m+half_eps,m,m),Point(zero,zero,zero)));
+	ASSERT_TRUE(lessEq(Point(m,m+half_eps,m),Point(zero,zero,zero)));
+	ASSERT_TRUE(lessEq(Point(m,m,m+half_eps),Point(zero,zero,zero)));
+	ASSERT_TRUE(lessEq(Point(m,m,m),Point(zero+half_eps,zero,zero)));
+	ASSERT_TRUE(lessEq(Point(m,m,m),Point(zero,zero+half_eps,zero)));
+	ASSERT_TRUE(lessEq(Point(m,m,m),Point(zero,zero,zero+half_eps)));
+
+	ASSERT_TRUE(lessEq(Point(zero+half_eps,zero,zero),Point(m,m,m)));
+	ASSERT_TRUE(lessEq(Point(zero,zero+half_eps,zero),Point(m,m,m)));
+	ASSERT_TRUE(lessEq(Point(zero,zero,zero+half_eps),Point(m,m,m)));
+	ASSERT_TRUE(lessEq(Point(zero,zero,zero),Point(m+half_eps,m,m)));
+	ASSERT_TRUE(lessEq(Point(zero,zero,zero),Point(m,m+half_eps,m)));
+	ASSERT_TRUE(lessEq(Point(zero,zero,zero),Point(m,m,m+half_eps)));
+
+	ASSERT_TRUE(lessEq(Point(half_eps+half_eps,half_eps,half_eps),Point(half_eps,half_eps,half_eps)));
+	ASSERT_TRUE(lessEq(Point(half_eps,half_eps+half_eps,zero),Point(half_eps,half_eps,half_eps)));
+	ASSERT_TRUE(lessEq(Point(zero,zero,zero+half_eps),Point(half_eps,half_eps,half_eps)));
+	ASSERT_TRUE(lessEq(Point(half_eps,half_eps,half_eps),Point(half_eps+half_eps,half_eps,half_eps)));
+	ASSERT_TRUE(lessEq(Point(half_eps,half_eps,half_eps),Point(half_eps,half_eps+half_eps,half_eps)));
+	ASSERT_TRUE(lessEq(Point(half_eps,half_eps,half_eps),Point(half_eps,half_eps,half_eps+half_eps)));
+
+	ASSERT_TRUE(lessEq(Point(10.0+half_eps,10.0,10.0),Point(10.0,10.0,10.0)));
+	ASSERT_TRUE(lessEq(Point(10.0,10.0+half_eps,10.0),Point(10.0,10.0,10.0)));
+	ASSERT_TRUE(lessEq(Point(10.0,10.0,10.0+half_eps),Point(10.0,10.0,10.0)));
+	ASSERT_TRUE(lessEq(Point(10.0,10.0,10.0),Point(10.0+half_eps,10.0,10.0)));
+	ASSERT_TRUE(lessEq(Point(10.0,10.0,10.0),Point(10.0,10.0+half_eps,10.0)));
+	ASSERT_TRUE(lessEq(Point(10.0,10.0,10.0),Point(10.0,10.0,10.0+half_eps)));
+}
+
+TEST(GeoLib, PointComparisonOperatorLessEq)
+{
+	const double my_eps(std::numeric_limits<double>::epsilon());
+	ASSERT_FALSE(Point(1.0+my_eps,1.0,1.0) <= Point(1.0,1.0,1.0));
+	ASSERT_FALSE(Point(1.0,1.0+my_eps,1.0) <= Point(1.0,1.0,1.0));
+	ASSERT_FALSE(Point(1.0,1.0,1.0+my_eps) <= Point(1.0,1.0,1.0));
+	ASSERT_TRUE(Point(1.0,1.0,1.0) <= Point(1.0+my_eps,1.0,1.0));
+	ASSERT_TRUE(Point(1.0,1.0,1.0) <= Point(1.0,1.0+my_eps,1.0));
+	ASSERT_TRUE(Point(1.0,1.0,1.0) <= Point(1.0,1.0,1.0+my_eps));
+
+	ASSERT_TRUE(Point(1.0-my_eps,1.0,1.0) <= Point(1.0,1.0,1.0));
+	ASSERT_TRUE(Point(1.0,1.0-my_eps,1.0) <= Point(1.0,1.0,1.0));
+	ASSERT_TRUE(Point(1.0,1.0,1.0-my_eps) <= Point(1.0,1.0,1.0));
+	ASSERT_FALSE(Point(1.0,1.0,1.0) <= Point(1.0-my_eps,1.0,1.0));
+	ASSERT_FALSE(Point(1.0,1.0,1.0) <= Point(1.0,1.0-my_eps,1.0));
+	ASSERT_FALSE(Point(1.0,1.0,1.0) <= Point(1.0,1.0,1.0-my_eps));
+
+	std::size_t n(10000);
+	srand ( time(NULL) );
+	for (std::size_t k(0); k<n; ++k) {
+		double random_val_x(((double)(rand()) / RAND_MAX - 0.5)); //real_dist(rng));
+		double random_val_y(((double)(rand()) / RAND_MAX - 0.5)); //real_dist(rng));
+		double random_val_z(((double)(rand()) / RAND_MAX - 0.5)); //real_dist(rng));
+
+		double big_x(random_val_x * std::numeric_limits<double>::max());
+		double big_y(random_val_y * std::numeric_limits<double>::max());
+		double big_z(random_val_z * std::numeric_limits<double>::max());
+
+		ASSERT_TRUE(Point(big_x-my_eps,big_y,big_z) <= Point(big_x,big_y,big_z));
+		ASSERT_TRUE(Point(big_x,big_y-my_eps,big_z) <= Point(big_x,big_y,big_z));
+		ASSERT_TRUE(Point(big_x,big_y,big_z-my_eps) <= Point(big_x,big_y,big_z));
+
+		ASSERT_TRUE(Point(big_x-my_eps,big_y-my_eps,big_z) <= Point(big_x,big_y,big_z));
+		ASSERT_TRUE(Point(big_x-my_eps,big_y,big_z-my_eps) <= Point(big_x,big_y,big_z));
+		ASSERT_TRUE(Point(big_x,big_y-my_eps,big_z-my_eps) <= Point(big_x,big_y,big_z));
+
+		ASSERT_TRUE(Point(big_x-my_eps,big_y-my_eps,big_z-my_eps) <= Point(big_x,big_y,big_z));
+
+		double small_x(random_val_x * std::numeric_limits<double>::epsilon());
+		double small_y(random_val_y * std::numeric_limits<double>::epsilon());
+		double small_z(random_val_z * std::numeric_limits<double>::epsilon());
+
+		ASSERT_TRUE(Point(small_x-my_eps,small_y,small_z) <= Point(small_x,small_y,small_z));
+		ASSERT_TRUE(Point(small_x,small_y-my_eps,small_z) <= Point(small_x,small_y,small_z));
+		ASSERT_TRUE(Point(small_x,small_y,small_z-my_eps) <= Point(small_x,small_y,small_z));
+
+		ASSERT_TRUE(Point(small_x-my_eps,small_y-my_eps,small_z) <= Point(small_x,small_y,small_z));
+		ASSERT_TRUE(Point(small_x-my_eps,small_y,small_z-my_eps) <= Point(small_x,small_y,small_z));
+		ASSERT_TRUE(Point(small_x,small_y-my_eps,small_z-my_eps) <= Point(small_x,small_y,small_z));
+
+		ASSERT_TRUE(Point(small_x-my_eps,small_y-my_eps,small_z-my_eps) <= Point(small_x,small_y,small_z));
+	}
+}