Merge pull request #13400 from Tytan:optimize_exposure_compensation

Optimize exposure compensation (#13400)

* Added perf test

* Optimized gains computation

* Use Eigen for gains calculation

modules/stitching/perf/perf_stich.cpp

@@ -13,6 +13,7 @@ using namespace perf;
 #define WORK_MEGAPIX 0.6
 
 typedef TestBaseWithParam<string> stitch;
+typedef TestBaseWithParam<int> stitchExposureCompensation;
 typedef TestBaseWithParam<tuple<string, string> > stitchDatasets;
 
 #ifdef HAVE_OPENCV_XFEATURES2D
@@ -20,6 +21,7 @@ typedef TestBaseWithParam<tuple<string, string> > stitchDatasets;
 #else
 #define TEST_DETECTORS testing::Values("orb", "akaze")
 #endif
+#define TEST_EXP_COMP_BS testing::Values(32, 16, 12, 10, 8)
 #define AFFINE_DATASETS testing::Values("s", "budapest", "newspaper", "prague")
 
 PERF_TEST_P(stitch, a123, TEST_DETECTORS)
@@ -58,6 +60,38 @@ PERF_TEST_P(stitch, a123, TEST_DETECTORS)
     SANITY_CHECK_NOTHING();
 }
 
+PERF_TEST_P(stitchExposureCompensation, a123, TEST_EXP_COMP_BS)
+{
+    Mat pano;
+
+    vector<Mat> imgs;
+    imgs.push_back( imread( getDataPath("stitching/a1.png") ) );
+    imgs.push_back( imread( getDataPath("stitching/a2.png") ) );
+    imgs.push_back( imread( getDataPath("stitching/a3.png") ) );
+
+    int bs = GetParam();
+
+    declare.time(30 * 10).iterations(10);
+
+    while(next())
+    {
+        Ptr<Stitcher> stitcher = Stitcher::create();
+        stitcher->setWarper(makePtr<SphericalWarper>());
+        stitcher->setRegistrationResol(WORK_MEGAPIX);
+        stitcher->setExposureCompensator(
+            makePtr<detail::BlocksGainCompensator>(bs, bs));
+
+        startTimer();
+        stitcher->stitch(imgs, pano);
+        stopTimer();
+    }
+
+    EXPECT_NEAR(pano.size().width, 1182, 50);
+    EXPECT_NEAR(pano.size().height, 682, 30);
+
+    SANITY_CHECK_NOTHING();
+}
+
 PERF_TEST_P(stitch, b12, TEST_DETECTORS)
 {
     Mat pano;

modules/stitching/src/exposure_compensate.cpp

@@ -41,6 +41,10 @@
 //M*/
 
 #include "precomp.hpp"
+#ifdef HAVE_EIGEN
+#include <Eigen/Core>
+#include <Eigen/Dense>
+#endif
 
 namespace cv {
 namespace detail {
@@ -80,6 +84,7 @@ void GainCompensator::feed(const std::vector<Point> &corners, const std::vector<
     const int num_images = static_cast<int>(images.size());
     Mat_<int> N(num_images, num_images); N.setTo(0);
     Mat_<double> I(num_images, num_images); I.setTo(0);
+    Mat_<bool> skip(num_images, 1); skip.setTo(true);
 
     //Rect dst_roi = resultRoi(corners, images);
     Mat subimg1, subimg2;
@@ -99,7 +104,19 @@ void GainCompensator::feed(const std::vector<Point> &corners, const std::vector<
                 submask2 = masks[j].first(Rect(roi.tl() - corners[j], roi.br() - corners[j])).getMat(ACCESS_READ);
                 intersect = (submask1 == masks[i].second) & (submask2 == masks[j].second);
 
-                N(i, j) = N(j, i) = std::max(1, countNonZero(intersect));
+                int intersect_count = countNonZero(intersect);
+                N(i, j) = N(j, i) = std::max(1, intersect_count);
+
+                // Don't compute Isums if subimages do not intersect anyway
+                if (intersect_count == 0)
+                    continue;
+
+                // Don't skip images that intersect with at least one other image
+                if (i != j)
+                {
+                    skip(i, 0) = false;
+                    skip(j, 0) = false;
+                }
 
                 double Isum1 = 0, Isum2 = 0;
                 for (int y = 0; y < roi.height; ++y)
@@ -123,22 +140,59 @@ void GainCompensator::feed(const std::vector<Point> &corners, const std::vector<
 
     double alpha = 0.01;
     double beta = 100;
+    int num_eq = num_images - countNonZero(skip);
 
-    Mat_<double> A(num_images, num_images); A.setTo(0);
-    Mat_<double> b(num_images, 1); b.setTo(0);
-    for (int i = 0; i < num_images; ++i)
+    Mat_<double> A(num_eq, num_eq); A.setTo(0);
+    Mat_<double> b(num_eq, 1); b.setTo(0);
+    for (int i = 0, ki = 0; i < num_images; ++i)
     {
-        for (int j = 0; j < num_images; ++j)
+        if (skip(i, 0))
+            continue;
+
+        for (int j = 0, kj = 0; j < num_images; ++j)
         {
-            b(i, 0) += beta * N(i, j);
-            A(i, i) += beta * N(i, j);
-            if (j == i) continue;
-            A(i, i) += 2 * alpha * I(i, j) * I(i, j) * N(i, j);
-            A(i, j) -= 2 * alpha * I(i, j) * I(j, i) * N(i, j);
+            if (skip(j, 0))
+                continue;
+
+            b(ki, 0) += beta * N(i, j);
+            A(ki, ki) += beta * N(i, j);
+            if (j != i)
+            {
+                A(ki, ki) += 2 * alpha * I(i, j) * I(i, j) * N(i, j);
+                A(ki, kj) -= 2 * alpha * I(i, j) * I(j, i) * N(i, j);
+            }
+            ++kj;
         }
+        ++ki;
     }
 
-    solve(A, b, gains_);
+    Mat_<double> l_gains;
+
+#ifdef HAVE_EIGEN
+    Eigen::MatrixXf eigen_A, eigen_b, eigen_x;
+    cv2eigen(A, eigen_A);
+    cv2eigen(b, eigen_b);
+
+    Eigen::LLT<Eigen::MatrixXf> solver(eigen_A);
+#if ENABLE_LOG
+    if (solver.info() != Eigen::ComputationInfo::Success)
+        LOGLN("Failed to solve exposure compensation system");
+#endif
+    eigen_x = solver.solve(eigen_b);
+
+    eigen2cv(eigen_x, l_gains);
+#else
+    solve(A, b, l_gains);
+#endif
+
+    gains_.create(num_images, 1);
+    for (int i = 0, j = 0; i < num_images; ++i)
+    {
+        if (skip(i, 0))
+            gains_.at<double>(i, 0) = 1;
+        else
+            gains_.at<double>(i, 0) = l_gains(j++, 0);
+    }
 
     LOGLN("Exposure compensation, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
 }