Add find_convex_hull using Graham scan

dlib/geometry/vector.h

@@ -1308,6 +1308,73 @@ 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::enable_if_t<std::is_same<T, point>::value || std::is_same<T, dpoint>::value, std::vector<T>>
+    find_convex_hull(std::vector<T>& points)
+    {
+        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 (p0.x() - a.x()) * (p0.x() - a.x()) + (p0.y() - a.y()) * (p0.y() - a.y()) <
+                           (p0.x() - b.x()) * (p0.x() - b.x()) + (p0.y() - b.y()) * (p0.y() - b.y());
+                default:
+                    throw std::logic_error("unreachable: invalid points_orientation");
+            }
+        });
+
+        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 <

dlib/geometry/vector_abstract.h

@@ -456,6 +456,21 @@ namespace dlib
               false if not.
     !*/
 
+// ----------------------------------------------------------------------------------------
+
+    template <typename T>
+    std::enable_if_t<std::is_same<T, point>::value || std::is_same<T, dpoint>::value, std::vector<T>>
+    find_convex_hull(
+        std::vector<T>& points
+    );
+    /*!
+        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 <

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;