Merge pull request #25607 from Fest1veNapkin:imgproc_approx_bounding_poly

Add a new function that approximates the polygon bounding a convex hull with a certain number of sides #25607

merge PR with <https://github.com/opencv/opencv_extra/pull/1179>

This PR is based on the paper [View Frustum Optimization To Maximize Object’s Image Area](https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=1fbd43f3827fffeb76641a9c5ab5b625eb5a75ba).

# Problem
I needed to reduce the number of vertices of the convex hull so that the additional area was minimal, andall vertices of the original contour enter the new contour.

![image](https://github.com/Fest1veNapkin/opencv/assets/98156294/efac35f6-b8f0-46ec-91e4-60800432620c)

![image](https://github.com/Fest1veNapkin/opencv/assets/98156294/2292d9d7-1c10-49c9-8489-23221b4b28f7)

# Description
Initially in the contour of n vertices, at each stage we consider the intersection points of the lines formed by each adjacent edges. Each of these intersection points will form a triangle with vertices through which lines pass. Let's choose a triangle with the minimum area and merge the two vertices at the intersection point. We continue until there are more vertices than the specified number of sides of the approximated polygon.
![image](https://github.com/Fest1veNapkin/opencv/assets/98156294/b87b21c4-112e-450d-a776-2a120048ca30)

# Complexity:
Using a std::priority_queue or std::set  time complexity is **(O(n\*ln(n))**, memory **O(n)**,
n - number of vertices in convex hull.

count of sides - the number of points by which we must reduce.
![image](https://github.com/Fest1veNapkin/opencv/assets/98156294/31ad5562-a67d-4e3c-bdc2-29f8b52caf88)

## Comment
If epsilon_percentage more 0, algorithm can return more values than _side_.
Algorithm returns OutputArray. If OutputArray.type() equals 0, algorithm returns values with InputArray.type().
New test uses image which are not in opencv_extra, needs to be added.

### Pull Request Readiness Checklist

See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request

- [ ] I agree to contribute to the project under Apache 2 License.
- [ ] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV
- [ ] The PR is proposed to the proper branch
- [ ] There is a reference to the original bug report and related work
- [ ] There is accuracy test, performance test and test data in opencv_extra repository, if applicable
      Patch to opencv_extra has the same branch name.
- [ ] The feature is well documented and sample code can be built with the project CMake

doc/opencv.bib

@@ -1467,3 +1467,12 @@
   volume = {60},
   journal = {ISPRS Journal of Photogrammetry and Remote Sensing}
 }
+@article{LowIlie2003,
+  author   = {Kok-Lim Low, Adrian Ilie},
+  year = {2003},
+  pages = {3-15},
+  title    = {View Frustum Optimization to Maximize Object's Image Area},
+  journal  = {Journal of Graphics, (GPU, & Game) Tools (JGT)},
+  volume = {8},
+  url = {https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=1fbd43f3827fffeb76641a9c5ab5b625eb5a75ba}
+}

modules/imgproc/include/opencv2/imgproc.hpp

@@ -4058,6 +4058,28 @@ CV_EXPORTS_W void approxPolyDP( InputArray curve,
                                 OutputArray approxCurve,
                                 double epsilon, bool closed );
 
+/** @brief Approximates a polygon with a convex hull with a specified accuracy and number of sides.
+
+The cv::approxPolyN function approximates a polygon with a convex hull
+so that the difference between the contour area of the original contour and the new polygon is minimal.
+It uses a greedy algorithm for contracting two vertices into one in such a way that the additional area is minimal.
+Straight lines formed by each edge of the convex contour are drawn and the areas of the resulting triangles are considered.
+Each vertex will lie either on the original contour or outside it.
+
+The algorithm based on the paper @cite LowIlie2003 .
+
+@param curve Input vector of a 2D points stored in std::vector or Mat, points must be float or integer.
+@param approxCurve Result of the approximation. The type is vector of a 2D point (Point2f or Point) in std::vector or Mat.
+@param nsides The parameter defines the number of sides of the result polygon.
+@param epsilon_percentage defines the percentage of the maximum of additional area.
+If it equals -1, it is not used. Otherwise algorighm stops if additional area is greater than contourArea(_curve) * percentage.
+If additional area exceeds the limit, algorithm returns as many vertices as there were at the moment the limit was exceeded.
+@param ensure_convex If it is true, algorithm creates a convex hull of input contour. Otherwise input vector should be convex.
+ */
+CV_EXPORTS_W void approxPolyN(InputArray curve, OutputArray approxCurve,
+                              int nsides, float epsilon_percentage = -1.0,
+                              bool ensure_convex = true);
+
 /** @brief Calculates a contour perimeter or a curve length.
 
 The function computes a curve length or a closed contour perimeter.

modules/imgproc/src/approx.cpp

@@ -39,6 +39,7 @@
 //
 //M*/
 #include "precomp.hpp"
+#include <queue>
 
 /****************************************************************************************\
 *                                  Chain Approximation                                   *
@@ -860,4 +861,240 @@ cvApproxPoly( const void* array, int header_size,
     return dst_seq;
 }
 
+enum class PointStatus : int8_t
+{
+    REMOVED = -1,
+    RECALCULATE = 0,
+    CALCULATED = 1
+};
+
+struct neighbours
+{
+    PointStatus pointStatus;
+    cv::Point2f point;
+    int next;
+    int prev;
+
+    explicit neighbours(int next_ = -1, int prev_ = -1, const cv::Point2f& point_ = { -1, -1 })
+    {
+        next = next_;
+        prev = prev_;
+        point = point_;
+        pointStatus = PointStatus::CALCULATED;
+    }
+};
+
+struct changes
+{
+    float area;
+    int vertex;
+    cv::Point2f intersection;
+
+    explicit changes(float area_, int vertex_, const cv::Point2f& intersection_)
+    {
+        area = area_;
+        vertex = vertex_;
+        intersection = intersection_;
+    }
+
+    bool operator < (const changes& elem) const
+    {
+        return (area < elem.area) || ((area == elem.area) && (vertex < elem.vertex));
+    }
+    bool operator > (const changes& elem) const
+    {
+        return (area > elem.area) || ((area == elem.area) && (vertex > elem.vertex));
+    }
+};
+
+/*
+  returns intersection point and extra area
+*/
+static void recalculation(std::vector<neighbours>& hull, int vertex_id, float& area_, float& x, float& y)
+{
+    cv::Point2f vertex = hull[vertex_id].point,
+        next_vertex = hull[hull[vertex_id].next].point,
+        extra_vertex_1 = hull[hull[vertex_id].prev].point,
+        extra_vertex_2 = hull[hull[hull[vertex_id].next].next].point;
+
+    cv::Point2f curr_edge = next_vertex - vertex,
+        prev_edge = vertex - extra_vertex_1,
+        next_edge = extra_vertex_2 - next_vertex;
+
+    float cross = prev_edge.x * next_edge.y - prev_edge.y * next_edge.x;
+    if (abs(cross) < 1e-8)
+    {
+        area_ = FLT_MAX;
+        x = -1;
+        y = -1;
+        return;
+    }
+
+    float t = (curr_edge.x * next_edge.y - curr_edge.y * next_edge.x) / cross;
+    cv::Point2f intersection = vertex + cv::Point2f(prev_edge.x * t, prev_edge.y * t);
+
+    float area = 0.5f * abs((next_vertex.x - vertex.x) * (intersection.y - vertex.y)
+        - (intersection.x - vertex.x) * (next_vertex.y - vertex.y));
+
+    area_ = area;
+    x = intersection.x;
+    y = intersection.y;
+}
+
+static void update(std::vector<neighbours>& hull, int vertex_id)
+{
+    neighbours& v1 = hull[vertex_id], & removed = hull[v1.next], & v2 = hull[removed.next];
+
+    removed.pointStatus = PointStatus::REMOVED;
+    v1.pointStatus = PointStatus::RECALCULATE;
+    v2.pointStatus = PointStatus::RECALCULATE;
+    hull[v1.prev].pointStatus = PointStatus::RECALCULATE;
+    v1.next = removed.next;
+    v2.prev = removed.prev;
+}
+
+/*
+    A greedy algorithm based on contraction of vertices for approximating a convex contour by a bounding polygon
+*/
+void cv::approxPolyN(InputArray _curve, OutputArray _approxCurve,
+    int nsides, float epsilon_percentage, bool ensure_convex)
+{
+    CV_INSTRUMENT_REGION();
+
+    CV_Assert(epsilon_percentage > 0 || epsilon_percentage == -1);
+    CV_Assert(nsides > 2);
+
+    if (_approxCurve.fixedType())
+    {
+        CV_Assert(_approxCurve.type() == CV_32FC2 || _approxCurve.type() == CV_32SC2);
+    }
+
+    Mat curve;
+    int depth = _curve.depth();
+
+    CV_Assert(depth == CV_32F || depth == CV_32S);
+
+    if (ensure_convex)
+    {
+        cv::convexHull(_curve, curve);
+    }
+    else
+    {
+        CV_Assert(isContourConvex(_curve));
+        curve = _curve.getMat();
+    }
+
+    CV_Assert((curve.cols == 1 && curve.rows >= nsides)
+        || (curve.rows == 1 && curve.cols >= nsides));
+
+    if (curve.rows == 1)
+    {
+        curve = curve.reshape(0, curve.cols);
+    }
+
+    std::vector<neighbours> hull(curve.rows);
+    int size = curve.rows;
+    std::priority_queue<changes, std::vector<changes>, std::greater<changes>> areas;
+    float extra_area = 0, max_extra_area = epsilon_percentage * static_cast<float>(contourArea(_curve));
+
+    if (curve.depth() == CV_32S)
+    {
+        for (int i = 0; i < size; ++i)
+        {
+            Point t = curve.at<cv::Point>(i, 0);
+            hull[i] = neighbours(i + 1, i - 1, Point2f(static_cast<float>(t.x), static_cast<float>(t.y)));
+        }
+    }
+    else
+    {
+        for (int i = 0; i < size; ++i)
+        {
+            Point2f t = curve.at<cv::Point2f>(i, 0);
+            hull[i] = neighbours(i + 1, i - 1, t);
+        }
+    }
+    hull[0].prev = size - 1;
+    hull[size - 1].next = 0;
+
+    if (size > nsides)
+    {
+        for (int vertex_id = 0; vertex_id < size; ++vertex_id)
+        {
+            float area, new_x, new_y;
+            recalculation(hull, vertex_id, area, new_x, new_y);
+
+            areas.push(changes(area, vertex_id, Point2f(new_x, new_y)));
+        }
+    }
+
+    while (size > nsides)
+    {
+        changes base = areas.top();
+        int vertex_id = base.vertex;
+
+        if (hull[vertex_id].pointStatus == PointStatus::REMOVED)
+        {
+            areas.pop();
+        }
+        else if (hull[vertex_id].pointStatus == PointStatus::RECALCULATE)
+        {
+            float area, new_x, new_y;
+            areas.pop();
+            recalculation(hull, vertex_id, area, new_x, new_y);
+
+            areas.push(changes(area, vertex_id, Point2f(new_x, new_y)));
+            hull[vertex_id].pointStatus = PointStatus::CALCULATED;
+        }
+        else
+        {
+            if (epsilon_percentage != -1)
+            {
+                extra_area += base.area;
+                if (extra_area > max_extra_area)
+                {
+                    break;
+                }
+            }
+
+            size--;
+            hull[vertex_id].point = base.intersection;
+            update(hull, vertex_id);
+        }
+    }
+
+    if (_approxCurve.fixedType())
+    {
+        depth = _approxCurve.depth();
+    }
+    _approxCurve.create(1, size, CV_MAKETYPE(depth, 2));
+    Mat buf = _approxCurve.getMat();
+    int last_free = 0;
+
+    if (depth == CV_32S)
+    {
+        for (int i = 0; i < curve.rows; ++i)
+        {
+            if (hull[i].pointStatus != PointStatus::REMOVED)
+            {
+                Point t = Point(static_cast<int>(round(hull[i].point.x)),
+                                static_cast<int>(round(hull[i].point.y)));
+
+                buf.at<Point>(0, last_free) = t;
+                last_free++;
+            }
+        }
+    }
+    else
+    {
+        for (int i = 0; i < curve.rows; ++i)
+        {
+            if (hull[i].pointStatus != PointStatus::REMOVED)
+            {
+                buf.at<Point2f>(0, last_free) = hull[i].point;
+                last_free++;
+            }
+        }
+    }
+}
+
 /* End of file. */

modules/imgproc/test/test_approxpoly.cpp

@@ -377,4 +377,79 @@ TEST(Imgproc_ApproxPoly, bad_epsilon)
     ASSERT_ANY_THROW(approxPolyDP(inputPoints, outputPoints, eps, false));
 }
 
+struct ApproxPolyN: public testing::Test
+{
+    void SetUp()
+    {
+        vector<vector<Point>> inputPoints = {
+            {  {87, 103}, {100, 112}, {96, 138}, {80, 169}, {60, 183}, {38, 176}, {41, 145}, {56, 118}, {76, 104} },
+            {  {196, 102}, {205, 118}, {174, 196}, {152, 207}, {102, 194}, {100, 175}, {131, 109} },
+            {  {372, 101}, {377, 119}, {337, 238}, {324, 248}, {240, 229}, {199, 214}, {232, 123}, {245, 103} },
+            {  {463, 86}, {563, 112}, {574, 135}, {596, 221}, {518, 298}, {412, 266}, {385, 164}, {462, 86} }
+        };
+
+        Mat image(600, 600, CV_8UC1, Scalar(0));
+
+        for (vector<Point>& polygon : inputPoints) {
+            polylines(image, { polygon }, true, Scalar(255), 1);
+        }
+
+        findContours(image, contours, RETR_LIST, CHAIN_APPROX_NONE);
+    }
+
+    vector<vector<Point>> contours;
+};
+
+TEST_F(ApproxPolyN, accuracyInt)
+{
+    vector<vector<Point>> rightCorners = {
+        { {72, 187}, {37, 176}, {42, 127}, {133, 64} },
+        { {168, 212}, {92, 192}, {131, 109}, {213, 100} },
+        { {72, 187}, {37, 176}, {42, 127}, {133, 64} },
+        { {384, 100}, {333, 251}, {197, 220}, {239, 103} },
+        { {168, 212}, {92, 192}, {131, 109}, {213, 100} },
+        { {333, 251}, {197, 220}, {239, 103}, {384, 100} },
+        { {542, 6}, {596, 221}, {518, 299}, {312, 236} },
+        { {596, 221}, {518, 299}, {312, 236}, {542, 6} }
+    };
+    EXPECT_EQ(rightCorners.size(), contours.size());
+
+    for (size_t i = 0; i < contours.size(); ++i) {
+        std::vector<Point> corners;
+        approxPolyN(contours[i], corners, 4, -1, true);
+        ASSERT_EQ(rightCorners[i], corners );
+    }
+}
+
+TEST_F(ApproxPolyN, accuracyFloat)
+{
+    vector<vector<Point2f>> rightCorners = {
+        { {72.f, 187.f}, {37.f, 176.f}, {42.f, 127.f}, {133.f, 64.f} },
+        { {168.f, 212.f}, {92.f, 192.f}, {131.f, 109.f}, {213.f, 100.f} },
+        { {72.f, 187.f}, {37.f, 176.f}, {42.f, 127.f}, {133.f, 64.f} },
+        { {384.f, 100.f}, {333.f, 251.f}, {197.f, 220.f}, {239.f, 103.f} },
+        { {168.f, 212.f}, {92.f, 192.f}, {131.f, 109.f}, {213.f, 100.f} },
+        { {333.f, 251.f}, {197.f, 220.f}, {239.f, 103.f}, {384.f, 100.f} },
+        { {542.f, 6.f}, {596.f, 221.f}, {518.f, 299.f}, {312.f, 236.f} },
+        { {596.f, 221.f}, {518.f, 299.f}, {312.f, 236.f}, {542.f, 6.f} }
+    };
+    EXPECT_EQ(rightCorners.size(), contours.size());
+
+    for (size_t i = 0; i < contours.size(); ++i) {
+        std::vector<Point2f> corners;
+        approxPolyN(contours[i], corners, 4, -1, true);
+        EXPECT_LT(cvtest::norm(rightCorners[i], corners, NORM_INF), .5f);
+    }
+}
+
+TEST_F(ApproxPolyN, bad_args)
+{
+    Mat contour(10, 1, CV_32FC2);
+    vector<vector<Point>> bad_contours;
+    vector<Point> corners;
+    ASSERT_ANY_THROW(approxPolyN(contour, corners, 0));
+    ASSERT_ANY_THROW(approxPolyN(contour, corners, 3, 0));
+    ASSERT_ANY_THROW(approxPolyN(bad_contours, corners, 4));
+}
+
 }} // namespace