Merge pull request #730 from alicevision/dev/imageMatchingMethods

Add Methods in imageMatching

src/aliceVision/geometry/Frustum.hpp

@@ -23,15 +23,17 @@ struct Frustum
 {
   Vec3 cones[5]; // camera centre and the 4 points that define the image plane
   Half_planes planes; // Define infinite frustum planes + 2 optional Near and Far Half Space
-  double z_near, z_far;
+  double z_near = -1., z_far = -1.;
   std::vector<Vec3> points;
 
-  Frustum() : z_near(-1.), z_far(-1.)  {}
+  Frustum() {}
 
   // Build a frustum from the image size, camera intrinsic and pose
-  Frustum(const int w, const int h, const Mat3 & K, const Mat3 & R, const Vec3 & C)
-    : z_near(-1.), z_far(-1.)
+  Frustum(const int w, const int h, const Mat3 & K, const Mat3 & R, const Vec3 & C, const double zNear = -1.0, const double zFar = -1.0)
+      : z_near(zNear)
+      , z_far(zFar)
   {
+    // supporting point are the points defined by the truncated cone
     const Mat3 Kinv = K.inverse();
     const Mat3 Rt = R.transpose();
 
@@ -48,40 +50,40 @@ struct Frustum
     planes.push_back( Half_plane_p(cones[0], cones[2], cones[3]) );
     planes.push_back( Half_plane_p(cones[0], cones[3], cones[4]) );
 
-    // supporting point for drawing is a normalized cone, since infinity cannot be represented
-    points = std::vector<Vec3>(&cones[0], &cones[0]+5);
-  }
-
-  Frustum(const int w, const int h, const Mat3 & K, const Mat3 & R, const Vec3 & C, const double zNear, const double zFar)
-  {
-    *this = Frustum(w,h,K,R,C);
-
-    // update near & far planes & clear set points
-    z_near = zNear;
-    z_far = zFar;
-    points.clear();
-    assert(zFar >= zNear);
-
     // Add Znear and ZFar half plane using the cam looking direction
     const Vec3 camLookDirection_n = R.row(2).normalized();
-    const double d_near = - zNear - camLookDirection_n.dot(C);
-    planes.push_back( Half_plane(camLookDirection_n, d_near) );
-
-    const double d_Far = zFar + camLookDirection_n.dot(C);
-    planes.push_back( Half_plane(-camLookDirection_n, d_Far) );
 
-    // supporting point are the points defined by the truncated cone
-    const Mat3 Kinv = K.inverse();
-    const Mat3 Rt = R.transpose();
-    points.push_back( Rt * (z_near * (Kinv * Vec3(0,0,1.0))) + C);
-    points.push_back( Rt * (z_near * (Kinv * Vec3(w,0,1.0))) + C);
-    points.push_back( Rt * (z_near * (Kinv * Vec3(w,h,1.0))) + C);
-    points.push_back( Rt * (z_near * (Kinv * Vec3(0,h,1.0))) + C);
-
-    points.push_back( Rt * (z_far * (Kinv * Vec3(0,0,1.0))) + C);
-    points.push_back( Rt * (z_far * (Kinv * Vec3(w,0,1.0))) + C);
-    points.push_back( Rt * (z_far * (Kinv * Vec3(w,h,1.0))) + C);
-    points.push_back( Rt * (z_far * (Kinv * Vec3(0,h,1.0))) + C);
+    if(zNear > 0)
+    {
+        const double d_near = - zNear - camLookDirection_n.dot(C);
+        planes.push_back( Half_plane(camLookDirection_n, d_near) );
+
+        points.push_back( Rt * (z_near * (Kinv * Vec3(0,0,1.0))) + C);
+        points.push_back( Rt * (z_near * (Kinv * Vec3(w,0,1.0))) + C);
+        points.push_back( Rt * (z_near * (Kinv * Vec3(w,h,1.0))) + C);
+        points.push_back( Rt * (z_near * (Kinv * Vec3(0,h,1.0))) + C);
+    }
+    else
+    {
+        points.push_back(cones[0]);
+    }
+    if(zFar > 0)
+    {
+        const double d_Far = zFar + camLookDirection_n.dot(C);
+        planes.push_back( Half_plane(-camLookDirection_n, d_Far) );
+
+        points.push_back( Rt * (z_far * (Kinv * Vec3(0,0,1.0))) + C);
+        points.push_back( Rt * (z_far * (Kinv * Vec3(w,0,1.0))) + C);
+        points.push_back( Rt * (z_far * (Kinv * Vec3(w,h,1.0))) + C);
+        points.push_back( Rt * (z_far * (Kinv * Vec3(0,h,1.0))) + C);
+    }
+    else
+    {
+        points.push_back(cones[1]);
+        points.push_back(cones[2]);
+        points.push_back(cones[3]);
+        points.push_back(cones[4]);
+    }
   }
 
   /// Test if two frustums intersect or not
@@ -107,6 +109,11 @@ struct Frustum
     return planes.size() == 6;
   }
 
+  bool isPartiallyTruncated() const
+  {
+    return planes.size() == 5;
+  }
+
   // Return the supporting frustum points (5 for the infinite, 8 for the truncated)
   const std::vector<Vec3> & frustum_points() const
   {

src/aliceVision/geometry/frustumIntersection_test.cpp

@@ -20,7 +20,7 @@
 using namespace aliceVision;
 using namespace aliceVision::geometry;
 using namespace aliceVision::geometry::halfPlane;
-using namespace std;
+
 
 //--
 // Camera frustum intersection unit test
@@ -81,6 +81,24 @@ BOOST_AUTO_TEST_CASE(intersection)
       for (int j = 0; j < iNviews; ++j)
         BOOST_CHECK(vec_frustum[i].intersect(vec_frustum[j]));
   }
+
+  //Test with partially truncated frustum
+  {
+    //Build frustum with near and far plane defined by a different value
+    std::vector<Frustum> vec_frustum;
+    for (int i=0; i < iNviews; ++i)
+    {
+      vec_frustum.push_back(
+        Frustum(principal_Point*2, principal_Point*2,
+                d._K[i], d._R[i], d._C[i], 0.1, -1.0));
+      BOOST_CHECK(vec_frustum[i].isPartiallyTruncated());
+    }
+
+    // Check that frustums have an overlap
+    for (int i = 0; i < iNviews; ++i)
+      for (int j = 0; j < iNviews; ++j)
+        BOOST_CHECK(vec_frustum[i].intersect(vec_frustum[j]));
+  }
 }
 
 BOOST_AUTO_TEST_CASE(empty_intersection)
@@ -169,3 +187,127 @@ BOOST_AUTO_TEST_CASE(empty_intersection)
     }
   }
 }
+
+void createPanoramaScene(double coeff, NViewDataSet& d, const int iNviews, const int iNbPoints, const int principal_Point)
+{
+    // Create a panorama scene
+    //--
+    // Create a panorama scene with a field of view
+    //  depending on a coefficient and the number of views
+    // 1 camera looks to iNviews different directions
+    //--
+
+    double field_of_view = (360.0*coeff)/iNviews;
+    double focalRatio = 0.5/(tan(0.5*degreeToRadian(field_of_view)));
+    double focalLengthPix=focalRatio*2*principal_Point;
+    NViewDatasetConfigurator config(focalLengthPix, focalLengthPix, principal_Point, principal_Point, 1, 0);
+    d._n = iNviews;
+    d._K.resize(iNviews);
+    d._R.resize(iNviews);
+    d._t.resize(iNviews);
+    d._C.resize(iNviews);
+
+    for (size_t i = 0; i < iNviews; ++i)
+    {
+      Vec3 camera_center(0., 0., 0.);
+      const double theta = i * 2 * M_PI / iNviews;
+      d._C[i] = camera_center;
+      // Circle
+      Vec3 lookdir(sin(theta), 0.0, cos(theta)); // Y axis UP
+
+      d._K[i] << config._fx,           0, config._cx,
+                          0,  config._fy, config._cy,
+                          0,           0,          1;
+      d._R[i] = LookAt(lookdir);  // Y axis UP
+      d._t[i] = -d._R[i] * camera_center; // [t]=[-RC] Cf HZ.
+    }
+
+}
+
+BOOST_AUTO_TEST_CASE(panorama_intersection)
+{
+  // Create partially truncated frustum
+  //--
+  // 1 camera looks to 4 different directions on a circle
+  // Create a panorama scene with a field of view
+  //  more than 90° for each view which means no overlap
+  //--
+
+  const int principal_Point = 500;
+  // Setup a panorama camera rig: cameras rotations around a nodal point
+  const int iNviews = 4;
+  const int iNbPoints = 6;
+  NViewDataSet d;
+  double coeff = 1.2; // overlap coefficient: more than 1 means overlap
+  // create panorama scene
+  createPanoramaScene(coeff, d, iNviews, iNbPoints, principal_Point);
+
+  //Test with partially truncated frustum
+  {
+    //Build frustum with near and far plane defined by a different value
+    std::vector<Frustum> vec_frustum;
+    for (int i=0; i < iNviews; ++i)
+    {
+      vec_frustum.push_back(
+        Frustum(principal_Point*2, principal_Point*2,
+                d._K[i], d._R[i], d._t[i], 0.1, -1.0));
+      BOOST_CHECK(vec_frustum[i].isPartiallyTruncated());
+    }
+
+    //Check that there is overlap between all frustums
+    for (int i = 0; i < iNviews; ++i)
+    {
+      int j = (i+1) % iNviews;
+      BOOST_CHECK(vec_frustum[i].intersect(vec_frustum[j]));
+
+      int k = (i-1+iNviews) % iNviews;
+      BOOST_CHECK(vec_frustum[i].intersect(vec_frustum[k]));
+    }
+  }
+}
+
+BOOST_AUTO_TEST_CASE(panorama_without_intersection)
+{
+  // Create partially truncated frustum
+  //--
+  // 1 camera looks to 4 different directions on a circle
+  // Create a panorama scene with a field of view
+  //  less than 90° for each view which means no overlap
+  //--
+
+  const int principal_Point = 500;
+  // Setup a panorama camera rig: cameras rotations around a nodal point
+  const int iNviews = 4;
+  const int iNbPoints = 6;
+  NViewDataSet d;
+  double coeff = 0.8;  // overlap coefficient: less than 1 means no overlap
+
+  //create panorama scene
+  createPanoramaScene(coeff, d, iNviews, iNbPoints, principal_Point);
+
+  //Test with partially truncated frustum
+  {
+    //Build frustum with near and far plane defined by a different value
+    std::vector<Frustum> vec_frustum;
+    for (int i=0; i < iNviews; ++i)
+    {
+      vec_frustum.push_back(
+        Frustum(principal_Point*2, principal_Point*2,
+                d._K[i], d._R[i], d._t[i], 0.1, -1.0));
+      BOOST_CHECK(vec_frustum[i].isPartiallyTruncated());
+    }
+
+    //Check that there is no overlap between all frustums
+    for (int i = 0; i < iNviews; ++i)
+    {
+      for (int j = 0; j < iNviews; ++j)
+      {
+        if(i == j)
+            continue;
+
+        BOOST_CHECK(!vec_frustum[i].intersect(vec_frustum[j]));
+      }
+    }
+  }
+}
+

src/aliceVision/matching/io.cpp

@@ -99,17 +99,22 @@ void filterMatchesByViews(
 
 void filterTopMatches(
   PairwiseMatches & allMatches,
-  const int limitNum)
+  const int limitNum,
+  const int minNum)
 {
-  if (limitNum <= 0)
+  if (limitNum <= 0 && minNum <=0)
     return;
+  if (limitNum > 0 && minNum > limitNum)
+    throw std::runtime_error("The minimum number of matches is higher than the maximum.");
 
   for(auto& matchesPerDesc: allMatches)
   {
     for(auto& matches: matchesPerDesc.second)
     {
       IndMatches& m = matches.second;
-      if (m.size() > limitNum)
+      if (minNum > 0 && m.size() < minNum)
+        m.clear();
+      else if (limitNum > 0 && m.size() > limitNum)
         m.erase(m.begin()+limitNum, m.end());
     }
   }
@@ -232,7 +237,8 @@ bool Load(
   const std::set<IndexT>& viewsKeysFilter,
   const std::vector<std::string>& folders,
   const std::vector<feature::EImageDescriberType>& descTypesFilter,
-  const int maxNbMatches)
+  const int maxNbMatches,
+  const int minNbMatches)
 {
   std::size_t nbLoadedMatchFiles = 0;
   const std::string pattern = "matches.txt";
@@ -271,8 +277,7 @@ bool Load(
   if(!descTypesFilter.empty())
     filterMatchesByDesc(matches, descTypesFilter);
 
-  if(maxNbMatches > 0)
-    filterTopMatches(matches, maxNbMatches);
+  filterTopMatches(matches, maxNbMatches, minNbMatches);
 
   ALICEVISION_LOG_TRACE("Matches per image pair (after filtering):");
   logMatches(matches);

src/aliceVision/matching/io.hpp

@@ -50,7 +50,8 @@ bool Load(PairwiseMatches& matches,
   const std::set<IndexT>& viewsKeysFilter,
   const std::vector<std::string>& folders,
   const std::vector<feature::EImageDescriberType>& descTypesFilter,
-  const int maxNbMatches = 0);
+  const int maxNbMatches = 0,
+  const int minNbMatches =0);
 
 /**
  * @brief Filter to keep only specific viewIds.

src/aliceVision/sfm/FrustumFilter.cpp

@@ -53,19 +53,9 @@ void FrustumFilter::initFrustum(const sfmData::SfMData& sfmData)
     if (cam == nullptr)
       continue;
 
-    if (!_bTruncated) // use infinite frustum
-    {
-      const Frustum f(
-        cam->w(), cam->h(), cam->K(),
-        pose.rotation(), pose.center());
-      frustum_perView[view->getViewId()] = f;
-    }
-    else // use truncated frustum with defined Near and Far planes
-    {
-      const Frustum f(cam->w(), cam->h(), cam->K(),
-        pose.rotation(), pose.center(), it->second.first, it->second.second);
-      frustum_perView[view->getViewId()] = f;
-    }
+    const Frustum f(cam->w(), cam->h(), cam->K(),
+      pose.rotation(), pose.center(), it->second.first, it->second.second);
+    frustum_perView[view->getViewId()] = f;
   }
 }
 

src/aliceVision/sfm/pipeline/ReconstructionEngine.cpp

@@ -4,8 +4,6 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file,
 // You can obtain one at https://mozilla.org/MPL/2.0/.
 
-#pragma once
-
 #include "ReconstructionEngine.hpp"
 
 #include <aliceVision/feature/RegionsPerView.hpp>

src/aliceVision/sfm/pipeline/pairwiseMatchesIO.hpp

@@ -35,6 +35,7 @@ inline bool loadPairwiseMatches(
     const std::vector<std::string>& folders,
     const std::vector<feature::EImageDescriberType>& descTypes,
     const int maxNbMatches = 0,
+    const int minNbMatches = 0,
     bool useOnlyMatchesFromFolder = false)
 {
   std::vector<std::string> matchesFolders;
@@ -47,7 +48,7 @@ inline bool loadPairwiseMatches(
   matchesFolders.insert(matchesFolders.end(), folders.begin(), folders.end());
 
   ALICEVISION_LOG_DEBUG("Loading matches");
-  if (!matching::Load(out_pairwiseMatches, sfmData.getViewsKeys(), matchesFolders, descTypes, maxNbMatches))
+  if (!matching::Load(out_pairwiseMatches, sfmData.getViewsKeys(), matchesFolders, descTypes, maxNbMatches, minNbMatches))
   {
     std::stringstream ss("Unable to read the matches file(s) from:\n");
     for(const std::string& folder : matchesFolders)

src/aliceVision/sfm/pipeline/sequential/ReconstructionEngine_sequentialSfM.cpp

@@ -236,11 +236,8 @@ std::size_t ReconstructionEngine_sequentialSfM::fuseMatchesIntoTracks()
     ALICEVISION_LOG_DEBUG("Track building");
     tracksBuilder.build(matches);
 
-    if(_params.useTrackFiltering)
-    {
-      ALICEVISION_LOG_DEBUG("Track filtering");
-      tracksBuilder.filter(_params.minInputTrackLength);
-    }
+    ALICEVISION_LOG_DEBUG("Track filtering");
+    tracksBuilder.filter(_params.filterTrackForks, _params.minInputTrackLength);
 
     ALICEVISION_LOG_DEBUG("Track export to internal structure");
     // build tracks with STL compliant type

src/aliceVision/sfm/pipeline/sequential/ReconstructionEngine_sequentialSfM.hpp

@@ -57,7 +57,7 @@ class ReconstructionEngine_sequentialSfM : public ReconstructionEngine
     double maxReprojectionError = 4.0;
     float minAngleInitialPair = 5.0f;
     float maxAngleInitialPair = 40.0f;
-    bool useTrackFiltering = true;
+    bool filterTrackForks = true;
     robustEstimation::ERobustEstimator localizerEstimator = robustEstimation::ERobustEstimator::ACRANSAC;
     double localizerEstimatorError = std::numeric_limits<double>::infinity();
     size_t localizerEstimatorMaxIterations = 4096;