Merge pull request #25469 from mshabunin:cpp-emd

imgproc: refactor EMD to reduce C-API usage #25469

- added more tests for EMD
- refactored to remove CvArr
- used BufferArea for memory allocations
- renamed functions and variables and formatted the code
- kept legacy functions intact in separate header

modules/imgproc/include/opencv2/imgproc/detail/legacy.hpp

@@ -23,6 +23,14 @@ CV_EXPORTS void findContours_legacy(InputArray image,
                                     int method,
                                     Point offset = Point());
 
+CV_EXPORTS float EMD_legacy( InputArray _signature1, InputArray _signature2,
+               int distType, InputArray _cost,
+               float* lowerBound, OutputArray _flow );
+
+CV_EXPORTS float wrapperEMD_legacy(InputArray _signature1, InputArray _signature2,
+               int distType, InputArray _cost,
+               Ptr<float> lowerBound, OutputArray _flow);
+
 #endif
 
 }  // namespace cv

modules/imgproc/src/emd.cpp

@@ -57,6 +57,7 @@
     ==========================================================================
 */
 #include "precomp.hpp"
+#include "opencv2/imgproc/detail/legacy.hpp"
 
 #define MAX_ITERATIONS 500
 #define CV_EMD_INF   ((float)1e20)
@@ -1147,7 +1148,7 @@ icvDistC( const float *x, const float *y, void *user_param )
 }
 
 
-float cv::EMD( InputArray _signature1, InputArray _signature2,
+float cv::EMD_legacy( InputArray _signature1, InputArray _signature2,
                int distType, InputArray _cost,
                float* lowerBound, OutputArray _flow )
 {
@@ -1171,7 +1172,7 @@ float cv::EMD( InputArray _signature1, InputArray _signature2,
                        _flow.needed() ? &_cflow : 0, lowerBound, 0 );
 }
 
-float cv::wrapperEMD(InputArray _signature1, InputArray _signature2,
+float cv::wrapperEMD_legacy(InputArray _signature1, InputArray _signature2,
                int distType, InputArray _cost,
                Ptr<float> lowerBound, OutputArray _flow)
 {

modules/imgproc/src/precomp.hpp

@@ -50,6 +50,8 @@
 #include "opencv2/core/private.hpp"
 #include "opencv2/core/ocl.hpp"
 #include "opencv2/core/hal/hal.hpp"
+#include "opencv2/core/check.hpp"
+#include "opencv2/core/utils/buffer_area.private.hpp"
 #include "opencv2/imgproc/hal/hal.hpp"
 #include "hal_replacement.hpp"
 

modules/imgproc/test/test_emd.cpp

@@ -1,93 +1,250 @@
-/*M///////////////////////////////////////////////////////////////////////////////////////
+// This file is part of OpenCV project.
-//
+// It is subject to the license terms in the LICENSE file found in the top-level directory
-//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+// of this distribution and at http://opencv.org/license.html
-//
-//  By downloading, copying, installing or using the software you agree to this license.
-//  If you do not agree to this license, do not download, install,
-//  copy or use the software.
-//
-//
-//                        Intel License Agreement
-//                For Open Source Computer Vision Library
-//
-// Copyright (C) 2000, Intel Corporation, all rights reserved.
-// Third party copyrights are property of their respective owners.
-//
-// Redistribution and use in source and binary forms, with or without modification,
-// are permitted provided that the following conditions are met:
-//
-//   * Redistribution's of source code must retain the above copyright notice,
-//     this list of conditions and the following disclaimer.
-//
-//   * Redistribution's in binary form must reproduce the above copyright notice,
-//     this list of conditions and the following disclaimer in the documentation
-//     and/or other materials provided with the distribution.
-//
-//   * The name of Intel Corporation may not be used to endorse or promote products
-//     derived from this software without specific prior written permission.
-//
-// This software is provided by the copyright holders and contributors "as is" and
-// any express or implied warranties, including, but not limited to, the implied
-// warranties of merchantability and fitness for a particular purpose are disclaimed.
-// In no event shall the Intel Corporation or contributors be liable for any direct,
-// indirect, incidental, special, exemplary, or consequential damages
-// (including, but not limited to, procurement of substitute goods or services;
-// loss of use, data, or profits; or business interruption) however caused
-// and on any theory of liability, whether in contract, strict liability,
-// or tort (including negligence or otherwise) arising in any way out of
-// the use of this software, even if advised of the possibility of such damage.
-//
-//M*/
 
+#include "opencv2/imgproc.hpp"
 #include "test_precomp.hpp"
 
+using namespace cv;
+using namespace std;
+
 namespace opencv_test { namespace {
 
-class CV_EMDTest : public cvtest::BaseTest
+//==============================================================================
+// Utility
+
+template <typename T>
+inline T sqr(T val)
 {
-public:
+    return val * val;
-    CV_EMDTest();
+}
-protected:
-    void run(int);
-};
 
+inline static float calcEMD(Mat w1, Mat w2, Mat& flow, int dist, int dims)
+{
+    float mass1 = 0.f, mass2 = 0.f, work = 0.f;
+    for (int i = 0; i < flow.rows; ++i)
+    {
+        mass1 += w1.at<float>(i, 0);
+        for (int j = 0; j < flow.cols; ++j)
+        {
+            if (i == 0)
+                mass2 += w2.at<float>(j, 0);
+            float dist_ = 0.f;
+            switch (dist)
+            {
+                case DIST_L1:
+                {
+                    for (int k = 1; k <= dims; ++k)
+                    {
+                        dist_ += abs(w1.at<float>(i, k) - w2.at<float>(j, k));
+                    }
+                    break;
+                }
+                case DIST_L2:
+                {
+                    for (int k = 1; k <= dims; ++k)
+                    {
+                        dist_ += sqr(w1.at<float>(i, k) - w2.at<float>(j, k));
+                    }
+                    dist_ = sqrt(dist_);
+                    break;
+                }
+                case DIST_C:
+                {
+                    for (int k = 1; k <= dims; ++k)
+                    {
+                        const float val = abs(w1.at<float>(i, k) - w2.at<float>(j, k));
+                        if (val > dist_)
+                            dist_ = val;
+                    }
+                    break;
+                }
+            }
+            const float weight = flow.at<float>(i, j);
+            work += dist_ * weight;
+        }
+    }
+    return work / max(mass1, mass2);
+}
+
+//==============================================================================
+
+TEST(Imgproc_EMD, regression)
+{
+    // input data
+    const float M = 10000;
+    Matx<float, 4, 1> w1 {50, 60, 50, 50};
+    Matx<float, 5, 1> w2 {30, 20, 70, 30, 60};
+    Matx<float, 4, 5> cost {16, 16, 13, 22, 17, 14, 14, 13, 19, 15,
+                            19, 19, 20, 23, M,  M,  0,  M,  0,  0};
+
+    // expected results
+    const double emd0 = 2460. / 210;
+    Matx<float, 4, 5> flow0 {0, 0, 50, 0, 0, 0, 0, 20, 0, 40, 30, 20, 0, 0, 0, 0, 0, 0, 30, 20};
 
-CV_EMDTest::CV_EMDTest()
+    // basic call with cost
     {
+        float emd = 0.f;
+        ASSERT_NO_THROW(emd = EMD(w1, w2, DIST_USER, cost));
+        EXPECT_NEAR(emd, emd0, 1e-6 * emd0);
     }
 
-void CV_EMDTest::run( int )
+    // basic call with cost and flow output
     {
-    int code = cvtest::TS::OK;
+        Mat flow;
-    const double success_error_level = 1e-6;
+        float emd = 0.f;
-    #define M 10000
+        ASSERT_NO_THROW(emd = EMD(w1, w2, DIST_USER, cost, nullptr, flow));
-    double emd0 = 2460./210;
+        EXPECT_NEAR(emd, emd0, 1e-6 * emd0);
-    static float cost[] =
+        EXPECT_MAT_NEAR(Mat(flow0), flow, 1e-6);
+    }
+    // no cost and DIST_USER - error
     {
-        16, 16, 13, 22, 17,
+        Mat flow;
-        14, 14, 13, 19, 15,
+        EXPECT_THROW(EMD(w1, w2, DIST_USER, noArray(), nullptr, flow), cv::Exception);
-        19, 19, 20, 23,  M,
+        EXPECT_THROW(EMD(w1, w2, DIST_USER), cv::Exception);
-        M ,  0,  M,  0,  0
+    }
-    };
+}
-    static float  w1[] = { 50, 60, 50, 50 },
-                  w2[] = { 30, 20, 70, 30, 60 };
-    Mat _w1(4, 1, CV_32F, w1);
-    Mat _w2(5, 1, CV_32F, w2);
-    Mat _cost(_w1.rows, _w2.rows, CV_32F, cost);
 
-    float emd = EMD( _w1, _w2, -1, _cost );
+TEST(Imgproc_EMD, distance_types)
-    if( fabs( emd - emd0 ) > success_error_level*emd0 )
 {
-        ts->printf( cvtest::TS::LOG,
+    // 1D (sum = 210)
-            "The computed distance is %.2f, while it should be %.2f\n", emd, emd0 );
+    Matx<float, 4, 2> w1 {50, 1, 60, 2, 50, 3, 50, 4};
-        code = cvtest::TS::FAIL_BAD_ACCURACY;
+    Matx<float, 5, 2> w2 {30, 1, 20, 2, 70, 3, 30, 4, 60, 5};
+
+    // 2D (sum = 210)
+    Matx<float, 4, 3> w3 {50, 0, 0, 60, 0, 1, 50, 1, 0, 50, 1, 1};
+    Matx<float, 5, 3> w4 {20, 0, 1, 70, 1, 0, 30, 1, 1, 60, 2, 2, 30, 3, 3};
+
+    // basic call with all distance types
+    {
+        const vector<DistanceTypes> good_types {DIST_L1, DIST_L2, DIST_C};
+        for (const auto& dt : good_types)
+        {
+            SCOPED_TRACE(cv::format("dt=%d", dt));
+            float emd = 0.f;
+            Mat flow;
+            // 1D
+            {
+                ASSERT_NO_THROW(emd = EMD(w1, w2, dt, noArray(), nullptr, flow));
+                const float emd0 = calcEMD(Mat(w1), Mat(w2), flow, dt, 1);
+                EXPECT_NEAR(emd0, emd, 1e-6);
+            }
+            // 2D
+            {
+                ASSERT_NO_THROW(emd = EMD(w3, w4, dt, noArray(), nullptr, flow));
+                const float emd0 = calcEMD(Mat(w3), Mat(w4), flow, dt, 2);
+                EXPECT_NEAR(emd0, emd, 1e-6);
+            }
+        }
     }
+}
+
+typedef testing::TestWithParam<int> Imgproc_EMD_dist;
+
+TEST_P(Imgproc_EMD_dist, random_flow_verify)
+{
+    const int dist = GetParam();
+    for (size_t iter = 0; iter < 100; ++iter)
+    {
+        SCOPED_TRACE(cv::format("iter=%zu", iter));
+        RNG& rng = TS::ptr()->get_rng();
+        const int dims = rng.uniform(1, 10);
+        Mat w1(rng.uniform(1, 10), dims + 1, CV_32FC1);
+        Mat w2(rng.uniform(1, 10), dims + 1, CV_32FC1);
+
+        // weights > 0
+        {
+            Mat w1_weights = w1.col(0);
+            Mat w2_weights = w2.col(0);
+            cvtest::randUni(rng, w1_weights, 0, 100);
+            cvtest::randUni(rng, w2_weights, 0, 100);
+        }
+
+        // coord
+        {
+            Mat w1_coord = w1.colRange(1, dims + 1);
+            Mat w2_coord = w2.colRange(1, dims + 1);
+            cvtest::randUni(rng, w1_coord, -10, +10);
+            cvtest::randUni(rng, w2_coord, -10, +10);
+        }
+
+        float emd1 = 0.f, emd2 = 0.f;
+        const float eps = 1e-5f;
+        Mat flow;
+        {
+            ASSERT_NO_THROW(emd1 = EMD(w1, w2, dist, noArray(), nullptr, flow));
+            const float emd0 = calcEMD(w1, w2, flow, dist, dims);
+            EXPECT_NEAR(emd0, emd1, eps);
+        }
+        {
+            ASSERT_NO_THROW(emd2 = EMD(w2, w1, dist, noArray(), nullptr, flow));
+            const float emd0 = calcEMD(w2, w1, flow, dist, dims);
+            EXPECT_NEAR(emd0, emd2, eps);
+        }
+        EXPECT_NEAR(emd1, emd2, eps);
+    }
+}
+
+INSTANTIATE_TEST_CASE_P(, Imgproc_EMD_dist, testing::Values(DIST_L1, DIST_L2, DIST_C));
 
-    if( code < 0 )
+
-        ts->set_failed_test_info( code );
+TEST(Imgproc_EMD, invalid)
+{
+    Matx<float, 4, 2> w1 {50, 1, 60, 2, 50, 3, 50, 4};
+    Matx<float, 5, 2> w2 {30, 1, 20, 2, 70, 3, 30, 4, 60, 5};
+
+    // empty signature
+    {
+        Mat empty;
+        EXPECT_THROW(EMD(empty, w2, DIST_USER), cv::Exception);
+        EXPECT_THROW(EMD(w1, empty, DIST_USER), cv::Exception);
+    }
+
+    // zero total weight, negative weight
+    {
+        Matx<float, 3, 1> wz {0, 0, 0};
+        Matx<float, 3, 2> wz1 {0, 1, 0, 2, 0, 3};
+        Matx<float, 3, 1> wn {0, 3, -2};
+        Matx<float, 3, 2> wn1 {0, 1, 3, 2, -2, 3};
+        EXPECT_THROW(EMD(wz, w2, DIST_USER), cv::Exception);
+        EXPECT_THROW(EMD(wz1, w2, DIST_USER), cv::Exception);
+        EXPECT_THROW(EMD(wn, w2, DIST_USER), cv::Exception);
+        EXPECT_THROW(EMD(wn1, w2, DIST_USER), cv::Exception);
+    }
+
+    // user distance type, but no cost matrix provided or is wrong
+    {
+        Mat cost(3, 3, CV_32FC1, Scalar::all(0)), cost8u(4, 5, CV_8UC1, Scalar::all(0)), empty;
+        EXPECT_THROW(EMD(w1, w2, DIST_USER, noArray()), cv::Exception);
+        EXPECT_THROW(EMD(w1, w2, DIST_USER, empty), cv::Exception);
+        EXPECT_THROW(EMD(w1, w2, DIST_USER, cost8u), cv::Exception);
+        EXPECT_THROW(EMD(w1, w2, DIST_USER, cost), cv::Exception);
     }
 
-TEST(Imgproc_EMD, regression) { CV_EMDTest test; test.safe_run(); }
+    // lower_bound is set together with cost
+    {
+        Mat cost(4, 5, CV_32FC1, Scalar::all(0));
+        float bound = 0.f;
+        EXPECT_THROW(EMD(w1, w2, DIST_USER, cost, &bound), cv::Exception);
+    }
+
+    // zero dimensions with non-user distance type
+    const vector<DistanceTypes> good_types {DIST_L1, DIST_L2, DIST_C};
+    for (const auto& dt : good_types)
+    {
+        SCOPED_TRACE(cv::format("dt=%d", dt));
+        Matx<float, 4, 1> w01 {20, 30, 40, 50};
+        Matx<float, 5, 1> w02 {20, 30, 40, 50, 10};
+        EXPECT_THROW(EMD(w01, w02, dt), cv::Exception);
+    }
+
+    // wrong distance type
+    const vector<DistanceTypes> bad_types {DIST_L12, DIST_FAIR, DIST_WELSCH, DIST_HUBER};
+    for (const auto& dt : bad_types)
+    {
+        SCOPED_TRACE(cv::format("dt=%d", dt));
+        EXPECT_THROW(EMD(w1, w2, dt), cv::Exception);
+    }
+}
 
-}} // namespace
+}}  // namespace opencv_test
-/* End of file. */

modules/imgproc/test/test_precomp.hpp

@@ -5,8 +5,11 @@
 #define __OPENCV_TEST_PRECOMP_HPP__
 
 #include "opencv2/ts.hpp"
+#include "opencv2/ts/ts_gtest.h"
+#include "opencv2/ts/ocl_test.hpp"
 #include "opencv2/imgproc.hpp"
 #include "opencv2/imgproc/imgproc_c.h"
+#include "opencv2/core.hpp"
 
 #include "opencv2/core/private.hpp"