fixed BA error calculation in opencv_stitching, added draft version of waveCorrect

modules/stitching/blenders.hpp

@@ -39,7 +39,7 @@ private:
 class MultiBandBlender : public Blender
 {
 public:
-    MultiBandBlender(int num_bands = 10) : num_bands_(num_bands) {}
+    MultiBandBlender(int num_bands = 7) : num_bands_(num_bands) {}
 
 private:
     cv::Point blend(const std::vector<cv::Mat> &src, const std::vector<cv::Point> &corners, const std::vector<cv::Mat> &masks,

modules/stitching/main.cpp

@@ -14,7 +14,8 @@ void printUsage()
     cout << "Rotation model images stitcher.\n" 
         << "Usage: opencv_stitching img1 img2 [...imgN]\n" 
         << "\t[--matchconf <0.0-1.0>]\n"
-        << "\t[--ba (ray_space|focal_ray_space)]\n"
+        << "\t[--ba (ray|focal_ray)]\n"
+        //<< "\t[--wavecorrect (no|yes)]\n"
         << "\t[--warp (plane|cylindrical|spherical)]\n" 
         << "\t[--seam (no|voronoi|graphcut)]\n" 
         << "\t[--blend (no|feather|multiband)]\n"
@@ -28,6 +29,7 @@ int main(int argc, char* argv[])
     vector<Mat> images;
     string result_name = "result.png";
     int ba_space = BundleAdjuster::RAY_SPACE;
+    bool wave_correct = false;
     int warp_type = Warper::SPHERICAL;
     bool user_match_conf = false;
     float match_conf = 0.55f;
@@ -55,9 +57,9 @@ int main(int argc, char* argv[])
         }
         else if (string(argv[i]) == "--ba")
         {
-            if (string(argv[i + 1]) == "ray_space")
+            if (string(argv[i + 1]) == "ray")
                 ba_space = BundleAdjuster::RAY_SPACE;
-            else if (string(argv[i + 1]) == "focal_ray_space")
+            else if (string(argv[i + 1]) == "focal_ray")
                 ba_space = BundleAdjuster::FOCAL_RAY_SPACE;
             else
             {
@@ -66,6 +68,19 @@ int main(int argc, char* argv[])
             }
             i++;
         }
+        //else if (string(argv[i]) == "--wavecorrect")
+        //{
+        //    if (string(argv[i + 1]) == "no")
+        //        wave_correct = false;
+        //    else if (string(argv[i + 1]) == "yes")
+        //        wave_correct = true;
+        //    else
+        //    {
+        //        cout << "Bad wave correct flag value\n";
+        //        return -1;
+        //    }
+        //    i++;
+        //}
         else if (string(argv[i]) == "--warp")
         {
             if (string(argv[i + 1]) == "plane")
@@ -164,6 +179,17 @@ int main(int argc, char* argv[])
     BundleAdjuster adjuster(ba_space);
     adjuster(images, features, pairwise_matches, cameras);
 
+    if (wave_correct)
+    {
+        LOGLN("Wave correcting...");
+        vector<Mat> rmats;
+        for (size_t i = 0; i < cameras.size(); ++i)
+            rmats.push_back(cameras[i].M);
+        waveCorrect(rmats);
+        for (size_t i = 0; i < cameras.size(); ++i)
+            cameras[i].M = rmats[i];
+    }
+
     // Find median focal length
     vector<double> focals;
     for (size_t i = 0; i < cameras.size(); ++i)

modules/stitching/motion_estimators.cpp

@@ -19,8 +19,8 @@ namespace
     public:
         inline CpuSurfFeaturesFinder() 
         {
-            detector_ = new SurfFeatureDetector(500.0);
-            extractor_ = new SurfDescriptorExtractor;
+            detector_ = new SurfFeatureDetector(500);
+            extractor_ = new SurfDescriptorExtractor();
         }
 
     protected:
@@ -491,17 +491,8 @@ void BundleAdjuster::estimate(const vector<Mat> &images, const vector<ImageFeatu
 
     edges_.clear();
     for (int i = 0; i < num_images_ - 1; ++i)
-    {
         for (int j = i + 1; j < num_images_; ++j)
-        {
-            int pair_idx = i * num_images_ + j;
-            const MatchesInfo& mi = pairwise_matches_[pair_idx];
-            float ni = static_cast<float>(mi.num_inliers);
-            float nf = static_cast<float>(mi.matches.size());
-            if (ni / (8.f + 0.3f * nf) > 1.f)
             edges_.push_back(make_pair(i, j));
-        }
-    }
 
     total_num_matches_ = 0;
     for (size_t i = 0; i < edges_.size(); ++i)
@@ -683,6 +674,43 @@ void BundleAdjuster::calcJacobian()
 }
 
 
+//////////////////////////////////////////////////////////////////////////////
+
+// TODO test on adobe/halfdome
+void waveCorrect(vector<Mat> &rmats)
+{
+    float data[9];
+    Mat r0(1, 3, CV_32F, data);
+    Mat r1(1, 3, CV_32F, data + 3);
+    Mat r2(1, 3, CV_32F, data + 6);
+    Mat R(3, 3, CV_32F, data);
+
+    Mat cov = Mat::zeros(3, 3, CV_32F);
+    for (size_t i = 0; i < rmats.size(); ++i)
+    {   
+        Mat r0 = rmats[i].col(0);
+        cov += r0 * r0.t();
+    }
+
+    SVD svd;
+    svd(cov, SVD::FULL_UV);
+    svd.vt.row(2).copyTo(r1);
+    if (determinant(svd.vt) < 0)
+        r1 *= -1;
+
+    Mat avgz = Mat::zeros(3, 1, CV_32F);
+    for (size_t i = 0; i < rmats.size(); ++i)
+        avgz += rmats[i].col(2);
+    avgz.t().cross(r1).copyTo(r0);
+    normalize(r0, r0);
+
+    r0.cross(r1).copyTo(r2);
+
+    for (size_t i = 0; i < rmats.size(); ++i)
+        rmats[i] = R * rmats[i];
+}
+
+
 //////////////////////////////////////////////////////////////////////////////
 
 void findMaxSpanningTree(int num_images, const vector<MatchesInfo> &pairwise_matches,

modules/stitching/motion_estimators.hpp

@@ -148,6 +148,9 @@ private:
 };
 
 
+void waveCorrect(std::vector<cv::Mat> &rmats);
+
+
 //////////////////////////////////////////////////////////////////////////////
 // Auxiliary functions