Merge pull request #1470 from IceRage:min_enclosing_triangle

modules/imgproc/doc/structural_analysis_and_shape_descriptors.rst

@@ -560,6 +560,41 @@ The function finds the four vertices of a rotated rectangle. This function is us
 
 
 
+minEnclosingTriangle
+----------------------
+Finds a triangle of minimum area enclosing a 2D point set and returns its area.
+
+.. ocv:function:: double minEnclosingTriangle( InputArray points, OutputArray triangle )
+
+.. ocv:pyfunction:: cv2.minEnclosingTriangle(points[, triangle]) -> retval, triangle
+
+    :param points: Input vector of 2D points with depth ``CV_32S`` or ``CV_32F``, stored in:
+
+            * ``std::vector<>`` or ``Mat`` (C++ interface)
+
+            * Nx2 numpy array (Python interface)
+
+    :param triangle: Output vector of three 2D points defining the vertices of the triangle. The depth of the OutputArray must be ``CV_32F``.
+
+The function finds a triangle of minimum area enclosing the given set of 2D points and returns its area. The output for a given 2D point set is shown in the image below. 2D points are depicted in *red* and the enclosing triangle in *yellow*.
+
+.. image:: pics/minenclosingtriangle.png
+    :height: 250px
+    :width: 250px
+    :alt: Sample output of the minimum enclosing triangle function
+
+The implementation of the algorithm is based on O'Rourke's [ORourke86]_ and Klee and Laskowski's [KleeLaskowski85]_ papers. O'Rourke provides a
+:math:`\theta(n)`
+algorithm for finding the minimal enclosing triangle of a 2D convex polygon with ``n`` vertices. Since the :ocv:func:`minEnclosingTriangle` function takes a 2D point set as input an additional preprocessing step of computing the convex hull of the 2D point set is required. The complexity of the :ocv:func:`convexHull` function is
+:math:`O(n log(n))` which is higher than
+:math:`\theta(n)`.
+Thus the overall complexity of the function is
+:math:`O(n log(n))`.
+
+.. note:: See ``opencv_source/samples/cpp/minarea.cpp`` for a usage example.
+
+
+
 minEnclosingCircle
 ----------------------
 Finds a circle of the minimum area enclosing a 2D point set.
@@ -672,6 +707,10 @@ See below a sample output of the function where each image pixel is tested again
 
 .. [Hu62] M. Hu. *Visual Pattern Recognition by Moment Invariants*, IRE Transactions on Information Theory, 8:2, pp. 179-187, 1962.
 
+.. [KleeLaskowski85] Klee, V. and Laskowski, M.C., *Finding the smallest triangles containing a given convex polygon*, Journal of Algorithms, vol. 6, no. 3, pp. 359-375 (1985)
+
+.. [ORourke86] O’Rourke, J., Aggarwal, A., Maddila, S., and Baldwin, M., *An optimal algorithm for finding minimal enclosing triangles*, Journal of Algorithms, vol. 7, no. 2, pp. 258-269 (1986)
+
 .. [Sklansky82] Sklansky, J., *Finding the Convex Hull of a Simple Polygon*. PRL 1 $number, pp 79-83 (1982)
 
 .. [Suzuki85] Suzuki, S. and Abe, K., *Topological Structural Analysis of Digitized Binary Images by Border Following*. CVGIP 30 1, pp 32-46 (1985)

modules/imgproc/include/opencv2/imgproc.hpp

@@ -1455,6 +1455,9 @@ CV_EXPORTS_W void boxPoints(RotatedRect box, OutputArray points);
 CV_EXPORTS_W void minEnclosingCircle( InputArray points,
                                       CV_OUT Point2f& center, CV_OUT float& radius );
 
+//! computes the minimal enclosing triangle for a set of points and returns its area
+CV_EXPORTS_W double minEnclosingTriangle( InputArray points, CV_OUT OutputArray triangle );
+
 //! matches two contours using one of the available algorithms
 CV_EXPORTS_W double matchShapes( InputArray contour1, InputArray contour2,
                                  int method, double parameter );

modules/imgproc/test/test_convhull.cpp

@@ -161,6 +161,22 @@ cvTsPointPolygonTest( CvPoint2D32f pt, const CvPoint2D32f* vv, int n, int* _idx=
     return result;
 }
 
+static cv::Point2f
+cvTsMiddlePoint(const cv::Point2f &a, const cv::Point2f &b)
+{
+    return cv::Point2f((a.x + b.x) / 2, (a.y + b.y) / 2);
+}
+
+static bool
+cvTsIsPointOnLineSegment(const cv::Point2f &x, const cv::Point2f &a, const cv::Point2f &b)
+{
+    double d1 = cvTsDist(CvPoint2D32f(x.x, x.y), CvPoint2D32f(a.x, a.y));
+    double d2 = cvTsDist(CvPoint2D32f(x.x, x.y), CvPoint2D32f(b.x, b.y));
+    double d3 = cvTsDist(CvPoint2D32f(a.x, a.y), CvPoint2D32f(b.x, b.y));
+
+    return (abs(d1 + d2 - d3) <= (1E-5));
+}
+
 
 /****************************************************************************************\
 *                              Base class for shape descriptor tests                     *
@@ -769,6 +785,145 @@ _exit_:
 }
 
 
+/****************************************************************************************\
+*                                   MinEnclosingTriangle Test                            *
+\****************************************************************************************/
+
+class CV_MinTriangleTest : public CV_BaseShapeDescrTest
+{
+public:
+    CV_MinTriangleTest();
+
+protected:
+    void run_func(void);
+    int validate_test_results( int test_case_idx );
+    std::vector<cv::Point2f> getTriangleMiddlePoints();
+
+    std::vector<cv::Point2f> convexPolygon;
+    std::vector<cv::Point2f> triangle;
+};
+
+
+CV_MinTriangleTest::CV_MinTriangleTest()
+{
+}
+
+std::vector<cv::Point2f> CV_MinTriangleTest::getTriangleMiddlePoints()
+{
+    std::vector<cv::Point2f> triangleMiddlePoints;
+
+    for (int i = 0; i < 3; i++) {
+        triangleMiddlePoints.push_back(cvTsMiddlePoint(triangle[i], triangle[(i + 1) % 3]));
+    }
+
+    return triangleMiddlePoints;
+}
+
+
+void CV_MinTriangleTest::run_func()
+{
+    std::vector<cv::Point2f> pointsAsVector;
+
+    cv::cvarrToMat(points).convertTo(pointsAsVector, CV_32F);
+
+    cv::minEnclosingTriangle(pointsAsVector, triangle);
+    cv::convexHull(pointsAsVector, convexPolygon, true, true);
+}
+
+
+int CV_MinTriangleTest::validate_test_results( int test_case_idx )
+{
+    bool errorEnclosed = false, errorMiddlePoints = false, errorFlush = true;
+    double eps = 1e-4;
+    int code = CV_BaseShapeDescrTest::validate_test_results( test_case_idx );
+
+#if 0
+    {
+    int n = 3;
+    double a = 8, c = 8, b = 100, d = 150;
+    CvPoint bp[4], *bpp = bp;
+    cvNamedWindow( "test", 1 );
+    IplImage* img = cvCreateImage( cvSize(500,500), 8, 3 );
+    cvZero(img);
+    for( i = 0; i < point_count; i++ )
+        cvCircle(img,cvPoint(cvRound(p[i].x*a+b),cvRound(p[i].y*c+d)), 3, CV_RGB(0,255,0), -1 );
+    for( i = 0; i < n; i++ )
+        bp[i] = cvPoint(cvRound(triangle[i].x*a+b),cvRound(triangle[i].y*c+d));
+    cvPolyLine( img, &bpp, &n, 1, 1, CV_RGB(255,255,0), 1, CV_AA, 0 );
+    cvShowImage( "test", img );
+    cvWaitKey();
+    cvReleaseImage(&img);
+    }
+#endif
+
+    int polygonVertices = (int) convexPolygon.size();
+
+    if (polygonVertices > 2) {
+        // Check if all points are enclosed by the triangle
+        for (int i = 0; (i < polygonVertices) && (!errorEnclosed); i++)
+        {
+            if (cv::pointPolygonTest(triangle, cv::Point2f(convexPolygon[i].x, convexPolygon[i].y), true) < (-eps))
+                errorEnclosed = true;
+        }
+
+        // Check if triangle edges middle points touch the polygon
+        std::vector<cv::Point2f> middlePoints = getTriangleMiddlePoints();
+
+        for (int i = 0; (i < 3) && (!errorMiddlePoints); i++)
+        {
+            bool isTouching = false;
+
+            for (int j = 0; (j < polygonVertices) && (!isTouching); j++)
+            {
+                if (cvTsIsPointOnLineSegment(middlePoints[i], convexPolygon[j],
+                                             convexPolygon[(j + 1) % polygonVertices]))
+                    isTouching = true;
+            }
+
+            errorMiddlePoints = (isTouching) ? false : true;
+        }
+
+        // Check if at least one of the edges is flush
+        for (int i = 0; (i < 3) && (errorFlush); i++)
+        {
+            for (int j = 0; (j < polygonVertices) && (errorFlush); j++)
+            {
+                if ((cvTsIsPointOnLineSegment(convexPolygon[j], triangle[i],
+                                              triangle[(i + 1) % 3])) &&
+                    (cvTsIsPointOnLineSegment(convexPolygon[(j + 1) % polygonVertices], triangle[i],
+                                              triangle[(i + 1) % 3])))
+                    errorFlush = false;
+            }
+        }
+
+        // Report any found errors
+        if (errorEnclosed)
+        {
+            ts->printf( cvtest::TS::LOG,
+            "All points should be enclosed by the triangle.\n" );
+            code = cvtest::TS::FAIL_BAD_ACCURACY;
+        }
+        else if (errorMiddlePoints)
+        {
+            ts->printf( cvtest::TS::LOG,
+            "All triangle edges middle points should touch the convex hull of the points.\n" );
+            code = cvtest::TS::FAIL_INVALID_OUTPUT;
+        }
+        else if (errorFlush)
+        {
+            ts->printf( cvtest::TS::LOG,
+            "At least one edge of the enclosing triangle should be flush with one edge of the polygon.\n" );
+            code = cvtest::TS::FAIL_INVALID_OUTPUT;
+        }
+    }
+
+    if ( code < 0 )
+        ts->set_failed_test_info( code );
+
+    return code;
+}
+
+
 /****************************************************************************************\
 *                                     MinEnclosingCircle Test                            *
 \****************************************************************************************/
@@ -1691,6 +1846,7 @@ void CV_PerimeterAreaSliceTest::run( int )
 
 TEST(Imgproc_ConvexHull, accuracy) { CV_ConvHullTest test; test.safe_run(); }
 TEST(Imgproc_MinAreaRect, accuracy) { CV_MinAreaRectTest test; test.safe_run(); }
+TEST(Imgproc_MinTriangle, accuracy) { CV_MinTriangleTest test; test.safe_run(); }
 TEST(Imgproc_MinCircle, accuracy) { CV_MinCircleTest test; test.safe_run(); }
 TEST(Imgproc_ContourPerimeter, accuracy) { CV_PerimeterTest test; test.safe_run(); }
 TEST(Imgproc_FitEllipse, accuracy) { CV_FitEllipseTest test; test.safe_run(); }

samples/cpp/minarea.cpp

@@ -8,12 +8,13 @@ using namespace std;
 
 static void help()
 {
-    cout << "This program demonstrates finding the minimum enclosing box or circle of a set\n"
+    cout << "This program demonstrates finding the minimum enclosing box, triangle or circle of a set\n"
-            "of points using functions: minAreaRect() minEnclosingCircle().\n"
+         << "of points using functions: minAreaRect() minEnclosingTriangle() minEnclosingCircle().\n"
-            "Random points are generated and then enclosed.\n"
+         << "Random points are generated and then enclosed.\n\n"
-            "Call:\n"
+         << "Press ESC, 'q' or 'Q' to exit and any other key to regenerate the set of points.\n\n"
-            "./minarea\n"
+         << "Call:\n"
-            "Using OpenCV v" << CV_VERSION << "\n" << endl;
+         << "./minarea\n"
+         << "Using OpenCV v" << CV_VERSION << "\n" << endl;
 }
 
 int main( int /*argc*/, char** /*argv*/ )
@@ -27,6 +28,8 @@ int main( int /*argc*/, char** /*argv*/ )
     {
         int i, count = rng.uniform(1, 101);
         vector<Point> points;
+
+        // Generate a random set of points
         for( i = 0; i < count; i++ )
         {
             Point pt;
@@ -36,23 +39,38 @@ int main( int /*argc*/, char** /*argv*/ )
             points.push_back(pt);
         }
 
+        // Find the minimum area enclosing bounding box
         RotatedRect box = minAreaRect(Mat(points));
 
+        // Find the minimum area enclosing triangle
+        vector<Point2f> triangle;
+
+        minEnclosingTriangle(points, triangle);
+
+        // Find the minimum area enclosing circle
         Point2f center, vtx[4];
         float radius = 0;
         minEnclosingCircle(Mat(points), center, radius);
         box.points(vtx);
 
         img = Scalar::all(0);
+
+        // Draw the points
         for( i = 0; i < count; i++ )
             circle( img, points[i], 3, Scalar(0, 0, 255), FILLED, LINE_AA );
 
+        // Draw the bounding box
         for( i = 0; i < 4; i++ )
             line(img, vtx[i], vtx[(i+1)%4], Scalar(0, 255, 0), 1, LINE_AA);
 
+        // Draw the triangle
+        for( i = 0; i < 3; i++ )
+            line(img, triangle[i], triangle[(i+1)%3], Scalar(255, 255, 0), 1, LINE_AA);
+
+        // Draw the circle
         circle(img, center, cvRound(radius), Scalar(0, 255, 255), 1, LINE_AA);
 
-        imshow( "rect & circle", img );
+        imshow( "Rectangle, triangle & circle", img );
 
         char key = (char)waitKey();
         if( key == 27 || key == 'q' || key == 'Q' ) // 'ESC'