From b20e97446b3e8cbb034d05efe5b37afe570e156a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Adri=C3=A0=20Arrufat?=
 <1671644+arrufat@users.noreply.github.com>
Date: Mon, 13 Nov 2023 10:10:47 +0900
Subject: [PATCH] Add find_convex_hull using Graham scan (#2889)

---
 dlib/geometry/vector.h          | 66 +++++++++++++++++++++++++++++++++
 dlib/geometry/vector_abstract.h | 16 ++++++++
 dlib/test/geometry.cpp          | 23 ++++++++++++
 3 files changed, 105 insertions(+)

diff --git a/dlib/geometry/vector.h b/dlib/geometry/vector.h
index 82af82b02..c5dbeeb8f 100644
--- a/dlib/geometry/vector.h
+++ b/dlib/geometry/vector.h
@@ -1308,6 +1308,72 @@ namespace dlib
         return (s0>0&&s1>0&&s2>0&&s3>0) || (s0<0&&s1<0&&s2<0&&s3<0);
     }
 
+// ----------------------------------------------------------------------------------------
+
+    namespace impl
+    {
+        enum class points_orientation
+        {
+            collinear,
+            clockwise,
+            couter_clockwise,
+        };
+
+        template <typename T>
+        std::enable_if_t<std::is_same<T, point>::value || std::is_same<T, dpoint>::value, points_orientation>
+        find_points_orientation(const T& a, const T& b, const T& c)
+        {
+            const auto v = a.x() * (b.y() - c.y()) + b.x() * (c.y() - a.y()) + c.x() * (a.y() - b.y());
+            if (v < 0) return points_orientation::clockwise;
+            if (v > 0) return points_orientation::couter_clockwise;
+            return points_orientation::collinear;
+        }
+    }
+
+    template <typename T>
+    std::vector<T> find_convex_hull(std::vector<T>& points)
+    {
+        static_assert(std::is_same<T, point>::value || std::is_same<T, dpoint>::value, "find_convex_hull() only works for 2D dlib::vector types.");
+        if (points.size() < 3)
+            return {};
+
+        // find the point with the lowest y coordinate, and the left-most in case of ties.
+        const auto p0 = *std::min_element(points.begin(), points.end(), [](const auto& a, const auto& b){
+            return std::make_pair(a.y(), a.x()) < std::make_pair(b.y(), b.x());
+        });
+
+        // sort the points by polar angle in clockwise order
+        std::sort(points.begin(), points.end(), [&p0](const auto& a, const auto& b){
+            switch (impl::find_points_orientation(p0, a, b))
+            {
+                case impl::points_orientation::clockwise:
+                    return true;
+                case impl::points_orientation::couter_clockwise:
+                    return false;
+                case impl::points_orientation::collinear:
+                    return length_squared(p0-a) < length_squared(p0-b);
+                default:
+                    DLIB_CASSERT(false, "this should be impossible");
+            }
+        });
+
+        std::vector<T> hull;
+        for (const auto& p : points)
+        {
+            while (hull.size() > 1 &&
+                impl::find_points_orientation(hull.at(hull.size() - 2), hull.back(), p) != impl::points_orientation::clockwise
+            )
+            {
+                hull.pop_back();
+            }
+            hull.push_back(p);
+        }
+        // If all the points were collinear, we'll have only two points in the hull, so we need to clear it.
+        if (hull.size() < 3)
+            hull.clear();
+        return hull;
+    }
+
 // ----------------------------------------------------------------------------------------
 
     template <
diff --git a/dlib/geometry/vector_abstract.h b/dlib/geometry/vector_abstract.h
index 6aa1b7ef6..5398a0c7e 100644
--- a/dlib/geometry/vector_abstract.h
+++ b/dlib/geometry/vector_abstract.h
@@ -456,6 +456,22 @@ namespace dlib
               false if not.
     !*/
 
+// ----------------------------------------------------------------------------------------
+
+    template <typename T>
+    std::vector<T> find_convex_hull(
+        std::vector<T>& points
+    );
+    /*!
+        requires
+            - T == dlib::point or dlib::dpoint
+        ensures
+            - If points.size() < 3: it returns an empty vector.
+            - Else: Finds the convex hull of points using the Graham scan algorithm.  That is,
+              the smallest convex shape that contains all points.  Moreover, in case all points
+              are collinear, that is, along the same line, it will also return an empty vector.
+    !*/
+
 // ----------------------------------------------------------------------------------------
 
     template <
diff --git a/dlib/test/geometry.cpp b/dlib/test/geometry.cpp
index 94a4aa6b5..dc04e656a 100644
--- a/dlib/test/geometry.cpp
+++ b/dlib/test/geometry.cpp
@@ -382,6 +382,28 @@ namespace
 
     }
 
+// ----------------------------------------------------------------------------------------
+
+    void test_find_convex_hull()
+    {
+        print_spinner();
+        std::vector<dpoint> points{
+            { 0.0, 0.0 },
+            { 1.0, 1.0 },
+            { 2.0, 2.0 },
+            { 3.0, 1.0 },
+            { 4.0, 0.0 },
+            { 2.0, 4.0 },
+            { 1.0, 3.0 },
+        };
+        const auto hull = find_convex_hull(points);
+        DLIB_TEST(hull.size() == 4);
+        DLIB_TEST(hull[0] == dpoint(0, 0));
+        DLIB_TEST(hull[1] == dpoint(1, 3));
+        DLIB_TEST(hull[2] == dpoint(2, 4));
+        DLIB_TEST(hull[3] == dpoint(4, 0));
+    }
+
 // ----------------------------------------------------------------------------------------
 
     void test_border_enumerator()
@@ -956,6 +978,7 @@ namespace
             test_find_similarity_transform2<float>(); 
             test_line();
             test_polygon();
+            test_find_convex_hull();
         }
     } a;