1 // This file is part of OpenMVG, an Open Multiple View Geometry C++ library.
2
3 // Copyright (c) 2015 Pierre MOULON.
4
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8
9 #include "openMVG/features/feature.hpp"
10 #include "openMVG/features/akaze/image_describer_akaze.hpp"
11 #include "openMVG/image/image_io.hpp"
12 #include "openMVG/matching/regions_matcher.hpp"
13 #include "openMVG/multiview/solver_homography_kernel.hpp"
14 #include "openMVG/robust_estimation/guided_matching.hpp"
15
16 #include "third_party/stlplus3/filesystemSimplified/file_system.hpp"
17 #include "third_party/cmdLine/cmdLine.h"
18
19 #include <string>
20 #include <iostream>
21
22 using namespace openMVG;
23 using namespace openMVG::image;
24 using namespace openMVG::matching;
25 using namespace std;
26
27 #include "openMVG/features/sift/SIFT_Anatomy_Image_Describer.hpp"
28 #include "nonFree/sift/SIFT_describer.hpp"
29
30 // Class to load images and ground truth homography matrices
31 // A reference image
32 // And a series of transformed images with the Homography mapping to the reference
33 class RepeatabilityDataset
34 {
35 public:
RepeatabilityDataset(const std::string & folderPath)36 explicit RepeatabilityDataset
37 (
38 const std::string& folderPath
39 ):
40 folderPath_(folderPath)
41 {
42 loadImages();
43 loadGroundTruthHs();
44 }
45
check() const46 bool check() const {
47 std::cout << "Dataset: " << folderPath_ << std::endl
48 << "#images: " << vec_image_.size() << "\n"
49 << "#homographies: " << vec_H_.size() << std::endl;
50 return !vec_H_.empty() && !vec_image_.empty() && vec_H_.size() == vec_image_.size();
51 }
52
image(size_t i) const53 const image::Image<RGBColor>& image(size_t i) const { return vec_image_[i]; }
H(size_t i) const54 const Mat3& H(size_t i) const { return vec_H_[i]; }
size() const55 size_t size() const { return vec_image_.size(); }
56
57 private:
58 /// Load the images of a directory
loadImages()59 bool loadImages()
60 {
61 std::cout << "Loading images of the dataset: " << stlplus::folder_part(folderPath_) << std::endl;
62 std::vector<std::string> vec_image_basename = stlplus::folder_wildcard( folderPath_, "*.ppm" );
63 if (vec_image_basename.empty())
64 vec_image_basename = stlplus::folder_wildcard( folderPath_, "*.pgm" );
65 if (vec_image_basename.empty())
66 return false;
67 sort(vec_image_basename.begin(), vec_image_basename.end());
68 vec_image_.resize(vec_image_basename.size());
69 for (size_t i = 0; i < vec_image_basename.size(); ++i)
70 {
71 const std::string path = stlplus::create_filespec(folderPath_, vec_image_basename[i]);
72 image::Image<RGBColor> imageRGB;
73 const bool bReadOk = image::ReadImage(path.c_str(), &imageRGB);
74 if ( bReadOk )
75 {
76 vec_image_[i] = imageRGB;
77 }
78 else
79 {
80 image::Image<unsigned char> imageGray;
81 if (image::ReadImage(path.c_str(), &imageGray))
82 {
83 image::ConvertPixelType(imageGray, &imageRGB);
84 vec_image_[i] = imageRGB;
85 }
86 else
87 {
88 std::cerr << "Error: unable to load image:\n" << path << std::endl;
89 return false;
90 }
91 }
92 }
93 return true;
94 }
95
96 /// Load the Homography related to each read images:
97 /// 0-> Identity,
98 /// 1-> H1to1p,
99 /// 2-> H1to2p, ...
loadGroundTruthHs()100 bool loadGroundTruthHs()
101 {
102 std::cout << "ground truth homographies of dataset: " << stlplus::folder_part(folderPath_) << std::endl;
103 vec_H_.resize(6);
104 for (int i = 0; i < 6; ++i)
105 {
106 if (i == 0)
107 {
108 vec_H_[i] = Mat3::Identity();
109 continue;
110 }
111
112 const std::string path = folderPath_ + "/H1to" + std::to_string(i+1) + "p";
113 std::ifstream f(path.c_str());
114 if (!f)
115 {
116 std::cerr << "Error: unable to load ground truth homography:\n"
117 << path << std::endl;
118 return false;
119 }
120 for (int k=0; k<9; ++k)
121 f >> vec_H_[i].data()[k];
122
123 vec_H_[i] /= vec_H_[i](2,2);
124
125 // Display
126 std::cout << "\n\n" << vec_H_[i] << std::endl;
127 }
128 return true;
129 }
130
131 private:
132 std::string folderPath_;
133
134 std::vector<image::Image<RGBColor>> vec_image_;
135 std::vector<Mat3> vec_H_;
136 };
137
138
139
140 /// Export point features based vector to matrices [(x,y)'T, (x,y)'T]
141 template< typename FeaturesT, typename MatT >
PointsToMat(const IndMatches & matches,const FeaturesT & vec_feats0,const FeaturesT & vec_feats1,MatT & m0,MatT & m1)142 void PointsToMat(
143 const IndMatches & matches,
144 const FeaturesT & vec_feats0,
145 const FeaturesT & vec_feats1,
146 MatT & m0,
147 MatT & m1)
148 {
149 using ValueT = typename FeaturesT::value_type; // Container type
150
151 m0.resize(2, matches.size());
152 m1.resize(2, matches.size());
153
154 for (size_t i = 0; i < matches.size(); ++i)
155 {
156 const ValueT & feat0 = vec_feats0[matches[i].i_];
157 m0.col(i) << feat0.x(), feat0.y();
158 const ValueT & feat1 = vec_feats1[matches[i].j_];
159 m1.col(i) << feat1.x(), feat1.y();
160 }
161 }
162
163 struct RepeatabilityResults_Matching
164 {
165 std::map<std::string, std::vector<double>> results;
166
exportToFileRepeatabilityResults_Matching167 bool exportToFile(const std::string & sFile, const std::string & sdatasetName) const
168 {
169 std::ofstream ofs(sFile, std::ofstream::out | std::ofstream::app);
170
171 if ( ! ofs.good() )
172 {
173 return false;
174 }
175
176 ofs << sdatasetName << "\n";
177 for ( const auto & val : results)
178 {
179 const std::string sParam = val.first;
180 const std::vector<double> & vec = val.second;
181 ofs << sParam << ";";
182 std::copy(vec.begin(), vec.end(), std::ostream_iterator<double>(ofs, ";"));
183 ofs << "\n";
184 }
185 ofs.close();
186
187 return true;
188 }
189 };
190
stringToEnum(const std::string & sPreset)191 features::EDESCRIBER_PRESET stringToEnum(const std::string & sPreset)
192 {
193 features::EDESCRIBER_PRESET preset;
194 if (sPreset == "NORMAL")
195 preset = features::NORMAL_PRESET;
196 else
197 if (sPreset == "HIGH")
198 preset = features::HIGH_PRESET;
199 else
200 if (sPreset == "ULTRA")
201 preset = features::ULTRA_PRESET;
202 else
203 preset = features::EDESCRIBER_PRESET(-1);
204 return preset;
205 }
206
207 //--
208 // Regions repeatability evaluation:
209 // - compare feature/descriptor matching repeatability on some dataset with known homography motions
210 // Must be run one of the dataset contained here:
211 // https://github.com/openMVG/Features_Repeatability
212 //
main(int argc,char ** argv)213 int main(int argc, char **argv)
214 {
215 CmdLine cmd;
216 //--
217 // Command line parameters
218 std::string sDataset_Path = "";
219 std::string sImage_Describer_Method = "SIFT";
220 std::string sFeature_Preset = "NORMAL";
221 bool bFeature_Repeatability = false;
222 bool bMatching_Repeatability = true;
223 cmd.add( make_option('i', sDataset_Path, "input_dataset") );
224 cmd.add( make_option('d', sImage_Describer_Method, "describer_method") );
225 cmd.add( make_option('p', sFeature_Preset, "describer_preset") );
226 cmd.add( make_option('f', bFeature_Repeatability, "feature_repeatability") );
227 cmd.add( make_option('m', bMatching_Repeatability, "matching_repeatability") );
228 //--
229
230 //--
231 // Command line parsing
232 try {
233 if (argc == 1) throw std::string("Invalid command line parameter.");
234 cmd.process(argc, argv);
235 } catch (const std::string& s) {
236 std::cerr << "Usage: " << argv[0] << '\n'
237 << "[-i|--input_dataset] the path to the datasets \n"
238 << "\n[Optional]\n"
239 << "[-d|--describer_method]\n"
240 << " (method to use to describe an image):\n"
241 << " SIFT (default),\n"
242 << " SIFT_ANATOMY,\n"
243 << " AKAZE_FLOAT: AKAZE with floating point descriptors.\n"
244 << "[-p|--describer_preset]\n"
245 << " (used to control the Image_describer configuration):\n"
246 << " NORMAL (default),\n"
247 << " HIGH,\n"
248 << " ULTRA: !!Can take long time!!\n"
249 << "[-f|--feature_repeatability]\n"
250 << " (used to control the Image_describer configuration):\n"
251 << " 0 (default),\n"
252 << " 1\n"
253 << "[-m|--matching_repeatability]\n"
254 << " 1 (default),\n"
255 << " 0\n"
256 << std::endl;
257
258 std::cerr << s << std::endl;
259 return EXIT_FAILURE;
260 }
261 //--
262
263 if (bFeature_Repeatability && bMatching_Repeatability)
264 {
265 std::cerr << "Only one repeatability test can be performed at a time." << std::endl;
266 return EXIT_FAILURE;
267 }
268
269 //--------------------------
270 // Evaluation parameters
271 //--------------------------
272 // - Precision radius to accept a point as a valid correspondence
273 const double m_dPrecision_robust = 2.5;
274 // - Nearest neighbor distance ratio to accept a photometric match between some descriptor
275 const double nndr = 0.8;
276
277 // List all subdirectories and for each one compute the repeatability
278 const std::vector<std::string> vec_dataset_path = stlplus::folder_subdirectories(sDataset_Path);
279 for (auto const& dataset_path : vec_dataset_path)
280 {
281 const std::string sPath = stlplus::create_filespec(sDataset_Path, dataset_path);
282 if (stlplus::is_file(sPath))
283 continue;
284
285 RepeatabilityDataset dataset(sPath);
286 if (dataset.check())
287 {
288 // 1. Compute regions
289 // 2. Test the repeatability (localization/overlap (accuracy))
290 // 3. Export data
291
292 using namespace openMVG::features;
293 std::unique_ptr<Image_describer> image_describer;
294 if (sImage_Describer_Method == "SIFT")
295 {
296 image_describer.reset(new SIFT_Image_describer);
297 }
298 else
299 if (sImage_Describer_Method == "SIFT_ANATOMY")
300 {
301 image_describer.reset(new SIFT_Anatomy_Image_describer);
302 }
303 else
304 if (sImage_Describer_Method == "AKAZE_FLOAT")
305 {
306 image_describer = AKAZE_Image_describer::create
307 (AKAZE_Image_describer::Params(AKAZE::Params(), AKAZE_MSURF));
308 }
309
310 if (!image_describer)
311 {
312 std::cerr << "Cannot create the designed Image_describer:"
313 << sImage_Describer_Method << "." << std::endl;
314 return EXIT_FAILURE;
315 }
316 else
317 {
318 if (!image_describer->Set_configuration_preset(stringToEnum(sFeature_Preset)))
319 {
320 std::cerr << "Preset configuration failed." << std::endl;
321 return EXIT_FAILURE;
322 }
323 }
324
325 // For each image computes the regions:
326 image::Image<unsigned char> imageGray;
327 Hash_Map<IndexT, std::unique_ptr<Regions>> map_regions;
328 for (size_t i = 0; i < dataset.size(); ++i)
329 {
330 image::ConvertPixelType(dataset.image(i), &imageGray);
331 image_describer->Describe(imageGray, map_regions[i]);
332 std::cout << "image: " << i << "\t #found features: " << map_regions[i]->RegionCount() << std::endl;
333 }
334
335 // Repeatability evaluation to the first image
336 // Evaluate the feature positions accuracy (descriptors are ignored)
337 if (bFeature_Repeatability)
338 {
339 RepeatabilityResults_Matching image_results;
340 for (size_t i = 1; i < dataset.size(); ++i)
341 {
342 if (map_regions.count(0) == 0 || map_regions.count(i) == 0)
343 continue;
344
345 const Regions * regions_0 = map_regions[0].get();
346 const Regions * regions_I = map_regions[i].get();
347 const PointFeatures pointsFeatures0 = regions_0->GetRegionsPositions();
348 const PointFeatures pointsFeaturesI = regions_I->GetRegionsPositions();
349
350 Mat x0, xI;
351 PointsToMat(pointsFeatures0, x0);
352 PointsToMat(pointsFeaturesI, xI);
353
354 IndMatches matches_0I;
355 geometry_aware::GuidedMatching
356 <Mat3, openMVG::homography::kernel::AsymmetricError>(
357 dataset.H(i).transpose(), x0, xI, Square(m_dPrecision_robust), matches_0I);
358
359 std::cout << "Feature repeatability Results" << "\n"
360 << "*******************************" << "\n"
361 << "# Keypoints 1: \t" << map_regions[0]->RegionCount() << "\n"
362 << "# Keypoints N: \t" << map_regions[i]->RegionCount() << "\n"
363 << "# Inliers: \t" << matches_0I.size() << "\n"
364 << "Inliers Ratio (%): \t" << matches_0I.size() / (float) map_regions[0]->RegionCount() << "\n"
365 << std::endl;
366
367 const std::vector<double> results = {
368 static_cast<double>( map_regions[0]->RegionCount() ) ,
369 static_cast<double>( map_regions[i]->RegionCount() ) ,
370 static_cast<double>( matches_0I.size() ) ,
371 static_cast<double>( matches_0I.size() / (float) map_regions[0]->RegionCount())
372 };
373 image_results.results[std::to_string(i)] = results;
374 }
375 image_results.exportToFile("repeatability_results.xls", stlplus::basename_part(sPath));
376 }
377
378 if (bMatching_Repeatability)
379 {
380 // Test the repeatability (matching (descriptor))
381 RepeatabilityResults_Matching image_results;
382 for (size_t i = 1; i < dataset.size(); ++i)
383 {
384 if (map_regions.count(0) == 0 || map_regions.count(i) == 0)
385 continue;
386
387 const Regions * regions_0 = map_regions[0].get();
388 const Regions * regions_I = map_regions[i].get();
389 const PointFeatures pointsFeatures0 = regions_0->GetRegionsPositions();
390 const PointFeatures pointsFeaturesI = regions_I->GetRegionsPositions();
391
392 matching::IndMatches putativesMatches;
393 matching::DistanceRatioMatch(
394 nndr, matching::BRUTE_FORCE_L2,
395 *regions_0, *regions_I,
396 putativesMatches);
397
398 Mat x0, xI;
399 PointsToMat(putativesMatches, pointsFeatures0, pointsFeaturesI, x0, xI);
400
401 IndMatches matches_0I;
402 geometry_aware::GuidedMatching
403 <Mat3, openMVG::homography::kernel::AsymmetricError>(
404 dataset.H(i).transpose(), x0, xI, Square(m_dPrecision_robust), matches_0I);
405
406 std::cout << "Feature matching repeatability Results" << "\n"
407 << "*******************************" << "\n"
408 << "** " << stlplus::basename_part(sPath) << " **" << "\n"
409 << "*******************************" << "\n"
410 << "# Keypoints 1: \t" << map_regions[0]->RegionCount() << "\n"
411 << "# Keypoints N: \t" << map_regions[i]->RegionCount() << "\n"
412 << "# Matches: \t" << putativesMatches.size() << "\n"
413 << "# Inliers: \t" << matches_0I.size() << "\n"
414 << "Inliers Ratio (%): \t" << matches_0I.size() / (float) putativesMatches.size() << "\n"
415 << std::endl;
416
417 const std::vector<double> results = {
418 static_cast<double>( map_regions[0]->RegionCount() ) ,
419 static_cast<double>( map_regions[i]->RegionCount() ) ,
420 static_cast<double>( putativesMatches.size() ) ,
421 static_cast<double>( matches_0I.size() ) ,
422 static_cast<double>( matches_0I.size() / (double) putativesMatches.size())
423 };
424 image_results.results[std::to_string(i)] = results;
425 }
426 image_results.exportToFile("repeatability_results.xls", stlplus::basename_part(sPath));
427 }
428 }
429 else
430 {
431 std::cerr << "Invalid dataset directory: " << dataset_path << std::endl;
432 return EXIT_FAILURE;
433 }
434 }
435 return EXIT_SUCCESS;
436 }
437