diff --git a/include/recti/interval.hpp b/include/recti/interval.hpp
index a9485b6..edc9a89 100644
--- a/include/recti/interval.hpp
+++ b/include/recti/interval.hpp
@@ -293,6 +293,18 @@ namespace recti {
             return x -= alpha;
         }
 
+        /**
+         * @brief Enlarge with
+         *
+         * @param[in] alpha
+         * @return interval<T>&
+         */
+        constexpr auto enlarge_with(const T& alpha) -> interval<T>& {
+            this->_lower -= alpha;
+            this->_upper += alpha;
+            return *this;
+        }
+
         ///@}
 
         /**
@@ -338,7 +350,9 @@ namespace recti {
          * @param[in] other
          * @return constexpr auto
          */
-        constexpr auto intersection(const T& other) const -> interval { return {other, other}; }
+        constexpr auto intersection_with(const T& other) const -> interval {
+            return {other, other};
+        }
 
         /**
          * @brief intersection with
@@ -346,7 +360,7 @@ namespace recti {
          * @param[in] other
          * @return constexpr auto
          */
-        constexpr auto intersection(const interval& other) const -> interval {
+        constexpr auto intersection_with(const interval& other) const -> interval {
             return {std::max(this->_lower, other._lower), std::min(this->_upper, other._upper)};
         }
 
@@ -419,10 +433,6 @@ namespace recti {
             *this = other = this->intersection(other);
             return 0;
         }
-
-        // ???
-        template <typename U1, typename U2>  //
-        friend inline constexpr auto min_dist_change(U1& lhs, U2& rhs);
     };
 #pragma pack(pop)
 
@@ -480,6 +490,15 @@ namespace recti {
         }
     }
 
+    template <typename U1, typename U2>  //
+    inline constexpr auto enlarge(U1 lhs, const U2& rhs) {
+        if constexpr (std::is_scalar_v<U1>) {
+            return interval<U1>{lhs - rhs, lhs + rhs};
+        } else {
+            return lhs.enlarge_with(rhs);
+        }
+    }
+
     template <typename U1, typename U2>  //
     inline constexpr auto min_dist_change(U1& lhs, U2& rhs) {
         if constexpr (std::is_scalar_v<U1>) {
diff --git a/include/recti/merge_obj.hpp b/include/recti/merge_obj.hpp
index 5022b1c..12cac41 100644
--- a/include/recti/merge_obj.hpp
+++ b/include/recti/merge_obj.hpp
@@ -3,6 +3,7 @@
 #include <algorithm>  // import std::min, std::max
 #include <cassert>
 #include <utility>  // import std::move
+
 #include "point.hpp"
 
 namespace recti {
@@ -23,7 +24,15 @@ namespace recti {
          * @param[in] x
          * @param[in] y
          */
-        constexpr merge_obj(const T1& x, const T2& y) :  point<T1, T2>{x + y, x - y} {}
+        constexpr merge_obj(T1&& x, T2&& y) noexcept : point<T1, T2>{std::move(x), std::move(y)} {}
+
+        // /**
+        //  * @brief Construct a new point object
+        //  *
+        //  * @param[in] x
+        //  * @param[in] y
+        //  */
+        // constexpr merge_obj(const T1& x, const T2& y) : point<T1, T2>{x + y, x - y} {}
 
         /**
          * @brief Add a vector (translation)
@@ -93,11 +102,36 @@ namespace recti {
             return std::max(min_dist(this->_x, other._x), min_dist(this->_y, other._y));
         }
 
+        template <typename R>  //
+        friend constexpr auto enlarge(const merge_obj& lhs, const R& alpha) {
+            auto x = enlarge(lhs.x(), alpha);
+            auto y = enlarge(lhs.y(), alpha);
+            return merge_obj<decltype(x), decltype(y)>{std::move(x), std::move(y)};
+        }
+
+        /**
+         * @brief overlap
+         *
+         * @tparam U1
+         * @tparam U2
+         * @param other
+         * @return true
+         * @return false
+         */
+        template <typename U1, typename U2>  //
+        [[nodiscard]] constexpr auto merge(const merge_obj<U1, U2>& other) const {
+            auto alpha = this->min_dist_with(other);
+            auto half = alpha / 2;
+            auto trr1 = enlarge(*this, half);
+            auto trr2 = ehlarge(other, alpha - half);
+            return trr1.intersection_with(trr2);
+        }
+
         // /**
         //  * @brief minimum distance with
-        //  * 
-        //  * @param[in] other 
-        //  * @return constexpr auto 
+        //  *
+        //  * @param[in] other
+        //  * @return constexpr auto
         //  */
         // [[nodiscard]] constexpr auto min_dist_change_with(merge_obj& other) {
         //     auto minDist = this->min_dist_with(other);
diff --git a/include/recti/point.hpp b/include/recti/point.hpp
index 6938a9a..29cf518 100644
--- a/include/recti/point.hpp
+++ b/include/recti/point.hpp
@@ -4,11 +4,19 @@
 #include <type_traits>  // import std::is_scalar_v
 #include <utility>      // import std::move
 
-#include "vector2.hpp"
 #include "interval.hpp"
+#include "vector2.hpp"
 
 namespace recti {
 
+    /**
+     * @brief Forward declaration
+     *
+     * @tparam U1
+     * @tparam U2
+     */
+    template <typename U1, typename U2> class point;
+
 #pragma pack(push, 1)
     /**
      * @brief 2D point
@@ -19,6 +27,14 @@ namespace recti {
     template <typename T1, typename T2 = T1> class point {
         using Self = point<T1, T2>;
 
+        /**
+         * @brief
+         *
+         * @tparam U1
+         * @tparam U2
+         */
+        template <typename U1, typename U2> friend class point;
+
       protected:
         T1 _x;  //!< x coordinate
         T2 _y;  //!< y coordinate
@@ -59,7 +75,7 @@ namespace recti {
          *
          * @return auto
          */
-        constexpr auto _tie() const { return std::tie(_x, _y); }
+        constexpr auto _tie() const { return std::tie(this->_x, this->_y); }
 
         /** @name Comparison operators
          *  definie ==, !=, <, >, <=, >=.
@@ -141,6 +157,8 @@ namespace recti {
         /**
          * @brief Greater than or equal to
          *
+         * @tparam U1
+         * @tparam U2
          * @param[in] x
          * @param[in] y
          * @return true
@@ -161,11 +179,13 @@ namespace recti {
         /**
          * @brief Add a vector (translation)
          *
-         * @tparam U
+         * @tparam U1
+         * @tparam U2
          * @param[in] rhs
          * @return Self&
          */
-        template <typename U> constexpr auto operator+=(const vector2<U>& rhs) -> Self& {
+        template <typename U1, typename U2> constexpr auto operator+=(const vector2<U1, U2>& rhs)
+            -> Self& {
             this->_x += rhs.x();
             this->_y += rhs.y();
             return *this;
@@ -174,11 +194,13 @@ namespace recti {
         /**
          * @brief Substract a vector (translation)
          *
-         * @tparam U
+         * @tparam U1
+         * @tparam U2
          * @param[in] rhs
          * @return Self&
          */
-        template <typename U> constexpr auto operator-=(const vector2<U>& rhs) -> Self& {
+        template <typename U1, typename U2> constexpr auto operator-=(const vector2<U1, U2>& rhs)
+            -> Self& {
             this->_x -= rhs.x();
             this->_y -= rhs.y();
             return *this;
@@ -187,26 +209,28 @@ namespace recti {
         /**
          * @brief Add
          *
-         * @tparam U
+         * @tparam U1
+         * @tparam U2
          * @param[in] x
          * @param[in] y
          * @return vector2<T>
          */
-        template <typename U>  //
-        friend constexpr auto operator+(point x, const vector2<U>& y) -> point {
+        template <typename U1, typename U2>  //
+        friend constexpr auto operator+(point x, const vector2<U1, U2>& y) -> point {
             return x += y;
         }
 
         /**
          * @brief Substract
          *
-         * @tparam U
+         * @tparam U1
+         * @tparam U2
          * @param[in] x
          * @param[in] y
          * @return vector2<T>
          */
-        template <typename U>  //
-        friend constexpr auto operator-(point x, const vector2<U>& y) -> point {
+        template <typename U1, typename U2>  //
+        friend constexpr auto operator-(point x, const vector2<U1, U2>& y) -> point {
             return x -= y;
         }
 
@@ -214,10 +238,12 @@ namespace recti {
          * @brief Different
          *
          * @param[in] rhs
-         * @return vector2
+         * @return constexpr auto
          */
-        constexpr auto operator-(const Self& rhs) const -> vector2<T1> {
-            return {this->x() - rhs.x(), this->y() - rhs.y()};
+        constexpr auto operator-(const Self& rhs) const {
+            auto x = this->x() - rhs.x();
+            auto y = this->y() - rhs.y();
+            return vector2<decltype(x), decltype(y)>{std::move(x), std::move(y)};
         }
 
         /**
@@ -280,7 +306,6 @@ namespace recti {
             return contain(this->x(), other.x()) && contain(this->y(), other.y());
         }
 
-
         /**
          * @brief overlap
          *
@@ -292,7 +317,7 @@ namespace recti {
          */
         template <typename U1, typename U2>  //
         [[nodiscard]] constexpr auto min_dist_with(const point<U1, U2>& other) const {
-            return min_dist(this->_x, other._x) + min_dist(this->_y, other._y);
+            return min_dist(this->x(), other.x()) + min_dist(this->y(), other.y());
         }
 
         /**
@@ -308,25 +333,31 @@ namespace recti {
         [[nodiscard]] constexpr auto min_dist_change_with(point<U1, U2>& other) {
             return min_dist_change(this->_x, other._x) + min_dist_change(this->_y, other._y);
         }
+
+        template <typename R>  //
+        friend constexpr auto enlarge(const point& lhs, const R& alpha) {
+            auto x = enlarge(lhs.x(), alpha);
+            auto y = enlarge(lhs.y(), alpha);
+            return point<decltype(x), decltype(y)>{std::move(x), std::move(y)};
+        }
+
+        /**
+         * @brief
+         *
+         * @tparam T1
+         * @tparam T2
+         * @tparam Stream
+         * @param[out] out
+         * @param[in] p
+         * @return Stream&
+         */
+        template <class Stream> friend auto operator<<(Stream& out, const point& p) -> Stream& {
+            out << "(" << p.x() << ", " << p.y() << ")";
+            return out;
+        }
     };
 #pragma pack(pop)
 
-    /**
-     * @brief
-     *
-     * @tparam T1
-     * @tparam T2
-     * @tparam Stream
-     * @param[out] out
-     * @param[in] p
-     * @return Stream&
-     */
-    template <typename T1, typename T2, class Stream>
-    auto operator<<(Stream& out, const point<T1, T2>& p) -> Stream& {
-        out << "(" << p.x() << ", " << p.y() << ")";
-        return out;
-    }
-
 #pragma pack(push, 1)
     /**
      * @brief 2D point
diff --git a/include/recti/vector2.hpp b/include/recti/vector2.hpp
index c636068..cf4a4c6 100644
--- a/include/recti/vector2.hpp
+++ b/include/recti/vector2.hpp
@@ -9,10 +9,10 @@ namespace recti {
      * @brief vector2
      *
      */
-    template <typename T = int> class vector2 {
+    template <typename T1, typename T2 = T1> class vector2 {
       private:
-        T _x;
-        T _y;
+        T1 _x;
+        T2 _y;
 
         auto _tie() const { return std::tie(_x, _y); }
 
@@ -23,7 +23,7 @@ namespace recti {
          * @param x
          * @param y
          */
-        constexpr vector2(T&& x, T&& y) noexcept : _x{std::move(x)}, _y{std::move(y)} {}
+        constexpr vector2(T1&& x, T2&& y) noexcept : _x{std::move(x)}, _y{std::move(y)} {}
 
         /**
          * @brief Construct a new vector2 object
@@ -31,38 +31,42 @@ namespace recti {
          * @param x
          * @param y
          */
-        constexpr vector2(const T& x, const T& y) : _x{x}, _y{y} {}
+        constexpr vector2(const T1& x, const T2& y) : _x{x}, _y{y} {}
 
         /**
          * @brief Construct a new vector2 object
          *
-         * @tparam U
+         * @tparam U1
+         * @tparam U2
          */
-        template <typename U> constexpr vector2(const vector2<U>& rhs) : _x(rhs.x()), _y(rhs.y()) {}
+        template <typename U1, typename U2> constexpr vector2(const vector2<U1, U2>& other)
+            : _x(other.x()), _y(other.y()) {}
 
         /**
          * @brief
          *
-         * @return constexpr const T&
+         * @return constexpr const T1&
          */
-        [[nodiscard]] constexpr auto x() const noexcept -> const T& { return this->_x; }
+        [[nodiscard]] constexpr auto x() const noexcept -> const T1& { return this->_x; }
 
         /**
          * @brief
          *
-         * @return constexpr const T&
+         * @return constexpr const T2&
          */
-        [[nodiscard]] constexpr auto y() const noexcept -> const T& { return this->_y; }
+        [[nodiscard]] constexpr auto y() const noexcept -> const T2& { return this->_y; }
 
         /**
          * @brief
          *
-         * @param[in] rhs
-         * @return T
+         * @tparam U1
+         * @tparam U2
+         * @param[in] other
+         * @return constexpr auto
          */
-        template <typename U>  //
-        [[nodiscard]] constexpr auto cross(const vector2<U>& rhs) const {
-            return this->_x * rhs._y - rhs._x * this->_y;
+        template <typename U1, typename U2>  //
+        [[nodiscard]] constexpr auto cross(const vector2<U1, U2>& other) const {
+            return this->_x * other._y - other._x * this->_y;
         }
 
         /** @name Comparison operators
@@ -73,79 +77,85 @@ namespace recti {
         /**
          * @brief Equal to
          *
-         * @tparam U
-         * @param[in] rhs
+         * @tparam U1
+         * @tparam U2
+         * @param[in] other
          * @return true
          * @return false
          */
-        template <typename U>  //
-        constexpr auto operator==(const vector2<U>& rhs) const -> bool {
-            return this->_tie() == rhs._tie();
+        template <typename U1, typename U2>  //
+        constexpr auto operator==(const vector2<U1, U2>& other) const -> bool {
+            return this->_tie() == other._tie();
         }
 
         /**
          * @brief Less than
          *
-         * @tparam U
-         * @param[in] rhs
+         * @tparam U1
+         * @tparam U2
+         * @param[in] other
          * @return true
          * @return false
          */
-        template <typename U>  //
-        constexpr auto operator<(const vector2<U>& rhs) const -> bool {
-            return this->_tie() < rhs._tie();
+        template <typename U1, typename U2>  //
+        constexpr auto operator<(const vector2<U1, U2>& other) const -> bool {
+            return this->_tie() < other._tie();
         }
 
         /**
          * @brief Not equal to
          *
-         * @tparam U
-         * @param[in] rhs
+         * @tparam U1
+         * @tparam U2
+         * @param[in] other
          * @return true
          * @return false
          */
-        template <typename U>  //
-        constexpr auto operator!=(const vector2<U>& rhs) const -> bool {
-            return !(*this == rhs);
+        template <typename U1, typename U2>  //
+        constexpr auto operator!=(const vector2<U1, U2>& other) const -> bool {
+            return !(*this == other);
         }
 
         /**
          * @brief Greater than
          *
-         * @tparam U
-         * @param[in] rhs
+         * @tparam U1
+         * @tparam U2
+         * @param[in] other
          * @return true
          * @return false
          */
-        template <typename U>  //
-        constexpr auto operator>(const vector2<U>& rhs) const -> bool {
-            return rhs < *this;
+        template <typename U1, typename U2>  //
+        constexpr auto operator>(const vector2<U1, U2>& other) const -> bool {
+            return other < *this;
         }
 
         /**
          * @brief Less than or euqual to
          *
-         * @tparam U
-         * @param[in] rhs
+         * @tparam U1
+         * @tparam U2
+         * @param[in] other
          * @return true
          * @return false
          */
-        template <typename U>  //
-        constexpr auto operator<=(const vector2<U>& rhs) const -> bool {
-            return !(rhs < *this);
+        template <typename U1, typename U2>  //
+        constexpr auto operator<=(const vector2<U1, U2>& other) const -> bool {
+            return !(other < *this);
         }
 
         /**
          * @brief Greater than or equal to
          *
-         * @tparam U
-         * @param[in] rhs
+         * @tparam U1
+         * @tparam U2
+         * @param[in] other
          * @return true
          * @return false
          */
-        template <typename U>  //
-        constexpr auto operator>=(const vector2<U>& rhs) const -> bool {
-            return !(*this < rhs);
+        template <typename U1, typename U2>  //
+        constexpr auto operator>=(const vector2<U1, U2>& other) const -> bool {
+            return !(*this < other);
         }
 
         ///@}
@@ -158,44 +168,48 @@ namespace recti {
         /**
          * @brief Negate
          *
-         * @return vector2<T>
+         * @return vector2
          */
-        constexpr auto operator-() const -> vector2<T> { return vector2<T>(-this->_x, -this->_y); }
+        constexpr auto operator-() const -> vector2 { return vector2(-this->_x, -this->_y); }
 
         /**
          * @brief Add
          *
-         * @tparam U
-         * @param[in] rhs
-         * @return vector2<T>&
+         * @tparam U1
+         * @tparam U2
+         * @param[in] other
+         * @return vector2&
          */
-        template <typename U> constexpr auto operator+=(const vector2<U>& rhs) -> vector2<T>& {
-            this->_x += rhs.x();
-            this->_y += rhs.y();
+        template <typename U1, typename U2> constexpr auto operator+=(const vector2<U1, U2>& other)
+            -> vector2& {
+            this->_x += other.x();
+            this->_y += other.y();
             return *this;
         }
 
         /**
          * @brief Substract
          *
-         * @tparam U
-         * @param[in] rhs
-         * @return vector2<T>&
+         * @tparam U1
+         * @tparam U2
+         * @param[in] other
+         * @return vector2&
          */
-        template <typename U>  //
-        constexpr auto operator-=(const vector2<U>& rhs) -> vector2<T>& {
-            this->_x -= rhs.x();
-            this->_y -= rhs.y();
+        template <typename U1, typename U2>  //
+        constexpr auto operator-=(const vector2<U1, U2>& other) -> vector2& {
+            this->_x -= other.x();
+            this->_y -= other.y();
             return *this;
         }
 
         /**
          * @brief Multiply
          *
+         * @tparam R
          * @param[in] alpha
-         * @return vector2<T>&
+         * @return vector2&
          */
-        constexpr auto operator*=(const T& alpha) -> vector2<T>& {
+        template <typename R> constexpr auto operator*=(const R& alpha) -> vector2& {
             this->_x *= alpha;
             this->_y *= alpha;
             return *this;
@@ -204,10 +218,11 @@ namespace recti {
         /**
          * @brief Divide
          *
+         * @tparam R
          * @param[in] alpha
-         * @return vector2<T>&
+         * @return vector2&
          */
-        constexpr auto operator/=(const T& alpha) -> vector2<T>& {
+        template <typename R> constexpr auto operator/=(const R& alpha) -> vector2& {
             this->_x /= alpha;
             this->_y /= alpha;
             return *this;
@@ -216,77 +231,83 @@ namespace recti {
         /**
          * @brief Add
          *
-         * @tparam U
+         * @tparam U1
+         * @tparam U2
          * @param[in] x
          * @param[in] y
-         * @return vector2<T>
+         * @return vector2
          */
-        template <typename U>  //
-        friend constexpr auto operator+(vector2<T> x, const vector2<U>& y) -> vector2<T> {
+        template <typename U1, typename U2>  //
+        friend constexpr auto operator+(vector2 x, const vector2<U1, U2>& y) -> vector2 {
             return x += y;
         }
 
         /**
          * @brief Substract
          *
-         * @tparam U
+         * @tparam U1
+         * @tparam U2
          * @param[in] x
          * @param[in] y
-         * @return vector2<T>
+         * @return vector2
          */
-        template <typename U>  //
-        friend constexpr auto operator-(vector2<T> x, const vector2<U>& y) -> vector2<T> {
+        template <typename U1, typename U2>  //
+        friend constexpr auto operator-(vector2 x, const vector2<U1, U2>& y) -> vector2 {
             return x -= y;
         }
 
         /**
          * @brief Multiply by a scalar
          *
+         * @tparam R
          * @param[in] x
          * @param[in] alpha
-         * @return vector2<T>
+         * @return vector2
          */
-        friend constexpr auto operator*(vector2<T> x, const T& alpha) -> vector2<T> {
+        template <typename R> friend constexpr auto operator*(vector2 x, const R& alpha)
+            -> vector2 {
             return x *= alpha;
         }
 
         /**
          * @brief Multiply (by a scalar)
          *
+         * @tparam R
          * @param[in] alpha
          * @param[in] x
-         * @return vector2<T>
+         * @return vector2
          */
-        friend constexpr auto operator*(const T& alpha, vector2<T> x) -> vector2<T> {
+        template <typename R> friend constexpr auto operator*(const R& alpha, vector2 x)
+            -> vector2 {
             return x *= alpha;
         }
 
         /**
          * @brief Divide (by a scalar)
          *
+         * @tparam R
          * @param[in] x
          * @param[in] alpha
-         * @return vector2<T>
+         * @return vector2
          */
-        friend constexpr auto operator/(vector2<T> x, const T& alpha) -> vector2<T> {
+        template <typename R> friend constexpr auto operator/(vector2 x, const R& alpha)
+            -> vector2 {
             return x /= alpha;
         }
 
         ///@}
-    };
 
-    /**
-     * @brief
-     *
-     * @tparam T
-     * @tparam Stream
-     * @param[out] out
-     * @param[in] v
-     * @return Stream&
-     */
-    template <typename T, class Stream> auto operator<<(Stream& out, const vector2<T>& v)
-        -> Stream& {
-        out << "{" << v.x() << ", " << v.y() << "}";
-        return out;
-    }
+        /**
+         * @brief
+         *
+         * @tparam Stream
+         * @param[out] out
+         * @param[in] v
+         * @return Stream&
+         */
+        template <class Stream> friend auto operator<<(Stream& out, const vector2& v) -> Stream& {
+            out << "{" << v.x() << ", " << v.y() << "}";
+            return out;
+        }
+    };
 }  // namespace recti
diff --git a/test/source/test_recti.cpp b/test/source/test_recti.cpp
index a57a1de..b3be353 100644
--- a/test/source/test_recti.cpp
+++ b/test/source/test_recti.cpp
@@ -52,9 +52,9 @@ TEST_CASE("Interval test") {
 
     CHECK(a.contains(4));
     CHECK(a.contains(8));
-    CHECK(a.intersection(8) == interval{8, 8});
+    CHECK(a.intersection_with(8) == interval{8, 8});
     CHECK(a.contains(b));
-    CHECK(a.intersection(b) == b);
+    CHECK(a.intersection_with(b) == b);
     CHECK(!b.contains(a));
     CHECK(a.overlaps(b));
     CHECK(b.overlaps(a));
@@ -134,8 +134,8 @@ TEST_CASE("Interval overlapping test") {
 }
 
 TEST_CASE("merge_obj test") {
-    auto r1 = merge_obj{4, 5};
-    auto r2 = merge_obj{7, 9};
+    auto r1 = merge_obj{4 + 5, 4 - 5};
+    auto r2 = merge_obj{7 + 9, 7 - 9};
     auto v = vector2{5, 6};
 
     CHECK(r1 != r2);
diff --git a/test/source/test_recti3d.cpp b/test/source/test_recti3d.cpp
new file mode 100644
index 0000000..82820bc
--- /dev/null
+++ b/test/source/test_recti3d.cpp
@@ -0,0 +1,151 @@
+#include <doctest/doctest.h>
+// #include <random>
+#include <cstdlib>
+#include <list>
+#include <ostream>
+#include <recti/halton_int.hpp>
+#include <recti/merge_obj.hpp>
+#include <recti/recti.hpp>
+#include <set>
+
+// using std::randint;
+using namespace recti;
+
+// template <typename T> struct my_point : point<T, T> { double data; };
+
+TEST_CASE("Point 3D test") {
+    auto a = point{point{40000, 80000}, 20000};
+    auto b = point{point{50000, 60000}, 10000};
+    auto v = (b - a) / 2;  // integer division
+
+    CHECK(a < b);
+    CHECK(a <= b);
+    CHECK(!(a == b));
+    CHECK(a != b);
+    CHECK(b > a);
+    CHECK(b >= a);
+    CHECK((a + v) + v == b);  // may not true due to integer division
+    CHECK((a - v) + v == a);
+
+    CHECK(a.flip_xy().flip_xy() == a);
+    // CHECK(a.flip_y().flip_y() == a);
+}
+
+// TEST_CASE("Interval test") {
+//     auto a = interval{4, 8};
+//     auto b = interval{5, 6};
+//     auto v = 3;
+
+//     CHECK(!(a < b));
+//     CHECK(!(b < a));
+//     CHECK(!(a > b));
+//     CHECK(!(b > a));
+//     CHECK(a <= b);
+//     CHECK(b <= a);
+//     CHECK(a >= b);
+//     CHECK(b >= a);
+
+//     CHECK(!(b == a));
+//     CHECK(b != a);
+
+//     CHECK((a - v) + v == a);
+
+//     CHECK(a.contains(4));
+//     CHECK(a.contains(8));
+//     CHECK(a.intersection(8) == interval{8, 8});
+//     CHECK(a.contains(b));
+//     CHECK(a.intersection(b) == b);
+//     CHECK(!b.contains(a));
+//     CHECK(a.overlaps(b));
+//     CHECK(b.overlaps(a));
+
+//     CHECK(min_dist(a, b) == 0);
+// }
+
+TEST_CASE("Rectangle 3D test") {
+    auto xrng1 = interval{40000, 80000};
+    auto yrng1 = interval{50000, 70000};
+    auto r1 = point{rectangle{xrng1, yrng1}, 1000};
+    auto xrng2 = interval{50000, 70000};
+    auto yrng2 = interval{60000, 60000};
+    auto r2 = point{rectangle{xrng2, yrng2}, 1000};
+    auto v = vector2{vector2{50000, 60000}, 0};
+    auto p1 = point{point{70000, 60000}, 1000};
+    auto p2 = point{point{70000, 60000}, 2000};
+
+    CHECK(r1 != r2);
+    CHECK((r1 - v) + v == r1);
+
+    // CHECK(r1 <= p);
+    CHECK(r1.contains(p1));
+    CHECK(!r1.contains(p2));
+    CHECK(r1.contains(r2));
+    CHECK(r1.overlaps(r2));
+    CHECK(overlap(r1, r2));
+
+    CHECK(r1.min_dist_with(r2) == 0);
+    CHECK(min_dist(r1, r2) == 0);
+
+    CHECK(r1.min_dist_with(p2) == p2.min_dist_with(r1));
+    // CHECK(min_dist(r1, p2) == min_dist(p2, r1));
+}
+
+// TEST_CASE("Segment test") {
+//     auto xrng1 = interval{4, 8};
+//     auto yrng1 = interval{5, 7};
+//     auto s1 = hsegment{xrng1, 6};
+//     auto s2 = vsegment{5, yrng1};
+
+//     CHECK(s1.overlaps(s2));
+// }
+
+// TEST_CASE("Interval overlapping test") {
+//     constexpr auto N = 20;
+//     auto lst = std::list<interval<unsigned int>>{};
+//     auto hgenX = vdcorput(3, 7);
+//     // auto hgenY = vdcorput(2, 11);
+
+//     for (auto i = 0; i != N; ++i) {
+//         for (auto j = 0; j != N; ++j) {
+//             auto x = hgenX();
+//             // auto y = hgenY();
+//             auto xrng = interval{x, x + 100};
+//             // auto yrng = interval{y, y + 100};
+//             // auto r = rectangle{xrng, yrng};
+//             // lst.push_back(r);
+//             lst.push_back(xrng);
+//         }
+//     }
+
+//     std::set<interval<unsigned int>> S;   // set of maximal non-overlapped rectangles
+//     std::list<interval<unsigned int>> L;  // list of the removed rectangles
+
+//     for (const auto& intvl : lst) {
+//         if (S.find(intvl) != S.end()) {
+//             L.push_back(intvl);
+//         } else {
+//             S.insert(intvl);
+//         }
+//     }
+
+//     // for (const auto& r : S)
+//     // {
+//     //     cout << "  \\draw " << r << ";\n";
+//     // }
+//     // for (const auto& r : L)
+//     // {
+//     //     cout << "  \\draw[color=red] " << r << ";\n";
+//     // }
+// }
+
+// TEST_CASE("merge_obj test") {
+//     auto r1 = merge_obj{4, 5};
+//     auto r2 = merge_obj{7, 9};
+//     auto v = vector2{5, 6};
+
+//     CHECK(r1 != r2);
+//     CHECK((r1 - v) + v == r1);
+//     CHECK(!overlap(r1, r2));
+//     CHECK(r1.min_dist_with(r2) == 7);
+//     CHECK(min_dist(r1, r2) == 7);
+// }