diff --git a/MathLib/GeometricBasics.h b/MathLib/GeometricBasics.h
index bc4b42fb5ce2ce1366825813b13217f6e95e62b5..f77ea72cb565f9cecf317f0dec0df08652ae8316 100644
--- a/MathLib/GeometricBasics.h
+++ b/MathLib/GeometricBasics.h
@@ -14,9 +14,7 @@
 
 namespace MathLib
 {
-template <typename T>
-class TemplatePoint;
-using Point3d = MathLib::TemplatePoint<double>;
+class Point3d;
 
 enum TriangleTest
 {
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index 1339e8ed3712c245348352eb7547d8a7b52e16fe..80a8377f56aff37baa39b6078d3c12879d5983d4 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -13,9 +13,7 @@
 
 namespace MathLib
 {
-template <typename T>
-class TemplatePoint;
-using Point3d = MathLib::TemplatePoint<double>;
+class Point3d;
 
 /**
  * calcProjPntToLineAndDists computes the orthogonal projection
diff --git a/MathLib/Point3d.cpp b/MathLib/Point3d.cpp
index 840d10ab0902ab237e63cdd52b2f8b1eb1453594..28b6f2394ddcc613eaeb960506595b70809adbba 100644
--- a/MathLib/Point3d.cpp
+++ b/MathLib/Point3d.cpp
@@ -12,5 +12,9 @@
 
 namespace MathLib
 {
+Point3d::Point3d() : x_({{0}}) {}
+
+Point3d::Point3d(std::array<double, 3> x) : x_(std::move(x)) {}
+
 extern const Point3d ORIGIN{{{0.0, 0.0, 0.0}}};
-}
+}  // namespace MathLib
diff --git a/MathLib/Point3d.h b/MathLib/Point3d.h
index 4bd0f9f55d0319814814cbd5c384ed8879214a0e..1375f8c82639337b840897b2eb1b85058679dd13 100644
--- a/MathLib/Point3d.h
+++ b/MathLib/Point3d.h
@@ -17,11 +17,117 @@
 
 #include "mathlib_export.h"
 
-#include "TemplatePoint.h"
-
 namespace MathLib
 {
-using Point3d = MathLib::TemplatePoint<double>;
+/**
+ * \ingroup MathLib
+ *
+ */
+class Point3d
+{
+public:
+    /** default constructor with zero coordinates */
+    Point3d();
+
+    /** constructor - constructs a Point3d object
+     *
+     * @param x std::array containing the coordinates of the point
+     */
+    explicit Point3d(std::array<double, 3> x);
+
+    /** virtual destructor */
+    virtual ~Point3d() = default;
+
+    Point3d(Point3d const&) = default;
+    Point3d& operator=(Point3d const&) = default;
+
+    /** \brief const access operator
+     *  The access to the point coordinates is like the access to a field. Code
+     * example: \code Point<double> point (1.0, 2.0, 3.0); double sqrNrm2 =
+     * point[0] * point[0] + point[1] * point[1] + point[2] + point[2]; \endcode
+     */
+    const double& operator[](std::size_t idx) const
+    {
+        assert(idx < 3);
+        return x_[idx];
+    }
+    /** \brief access operator (see book Effektiv C++ programmieren -
+     * subsection 1.3.2 ). \sa const T& operator[] (std::size_t idx) const
+     */
+    double& operator[](std::size_t idx)
+    {
+        return const_cast<double&>(static_cast<const Point3d&>(*this)[idx]);
+    }
+
+    /** returns an array containing the coordinates of the point */
+    const double* getCoords() const { return x_.data(); }
+
+    double* getCoords() { return x_.data(); }
+
+private:
+    std::array<double, 3> x_;
+};
+
+inline bool operator<(Point3d const& a, Point3d const& b)
+{
+    for (std::size_t i = 0; i < 3; ++i)
+    {
+        if (a[i] > b[i])
+        {
+            return false;
+        }
+        if (a[i] < b[i])
+        {
+            return true;
+        }
+
+        // continue with next dimension, because a[0] == b[0]
+    }
+
+    // The values in all dimensions are equal.
+    return false;
+}
+
+/**
+ * Lexicographic comparison of points taking an epsilon into account.
+ *
+ * @param a first input point.
+ * @param b second input point.
+ * @param eps tolerance used in comparison of coordinates.
+ *
+ * @return true, if a is smaller then or equal to b according to the following
+ * test \f$ |a_i - b_i| > \epsilon \cdot \min (|a_i|, |b_i|) \f$ \b and
+ * \f$  |a_i - b_i| > \epsilon \f$ for all coordinates \f$ 0 \le i < 3 \f$.
+ */
+bool inline lessEq(Point3d const& a, Point3d const& b,
+                   double eps = std::numeric_limits<double>::epsilon())
+{
+    auto coordinateIsLargerEps = [&eps](double const u, double const v) -> bool
+    {
+        return std::abs(u - v) > eps * std::min(std::abs(v), std::abs(u)) &&
+               std::abs(u - v) > eps;
+    };
+
+    for (std::size_t i = 0; i < 3; ++i)
+    {
+        // test a relative and an absolute criterion
+        if (coordinateIsLargerEps(a[i], b[i]))
+        {
+            return static_cast<bool>(a[i] <= b[i]);
+        }
+        // a[i] ~= b[i] up to an epsilon. Compare next dimension.
+    }
+
+    // all coordinates are equal up to an epsilon.
+    return true;
+}
+
+/** overload the output operator for class Point */
+inline std::ostream& operator<<(std::ostream& os, const Point3d& p)
+{
+    os << p[0] << " " << p[1] << " " << p[2];
+    return os;
+}
 
 extern MATHLIB_EXPORT const Point3d ORIGIN;
 /**
@@ -34,9 +140,11 @@ template <typename MATRIX>
 inline MathLib::Point3d operator*(MATRIX const& mat, MathLib::Point3d const& p)
 {
     MathLib::Point3d new_p;
-    for (std::size_t i(0); i<3; ++i) {
-        for (std::size_t j(0); j<3; ++j) {
-            new_p[i] += mat(i,j)*p[j];
+    for (std::size_t i(0); i < 3; ++i)
+    {
+        for (std::size_t j(0); j < 3; ++j)
+        {
+            new_p[i] += mat(i, j) * p[j];
         }
     }
     return new_p;
@@ -44,20 +152,28 @@ inline MathLib::Point3d operator*(MATRIX const& mat, MathLib::Point3d const& p)
 
 /** Computes the squared dist between the two points p0 and p1.
  */
-inline
-double sqrDist(MathLib::Point3d const& p0, MathLib::Point3d const& p1)
+inline double sqrDist(MathLib::Point3d const& p0, MathLib::Point3d const& p1)
 {
     return (p0[0] - p1[0]) * (p0[0] - p1[0]) +
            (p0[1] - p1[1]) * (p0[1] - p1[1]) +
            (p0[2] - p1[2]) * (p0[2] - p1[2]);
 }
 
+/** Equality of Point3d's up to an epsilon.
+ */
+inline bool operator==(Point3d const& a, Point3d const& b)
+{
+    auto const sqr_dist(sqrDist(a, b));
+    auto const eps = std::numeric_limits<double>::epsilon();
+    return (sqr_dist < eps * eps);
+}
+
 /// Computes the squared distance between the orthogonal projection of the two
 /// points \c p0 and \c p1 onto the \f$xy\f$-plane.
-inline
-double sqrDist2d(MathLib::Point3d const& p0, MathLib::Point3d const& p1)
+inline double sqrDist2d(MathLib::Point3d const& p0, MathLib::Point3d const& p1)
 {
-    return (p0[0]-p1[0])*(p0[0]-p1[0]) + (p0[1]-p1[1])*(p0[1]-p1[1]);
+    return (p0[0] - p1[0]) * (p0[0] - p1[0]) +
+           (p0[1] - p1[1]) * (p0[1] - p1[1]);
 }
 
-} // end namespace MathLib
+}  // end namespace MathLib