38 #include <visp3/core/vpConfig.h>
40 #if defined(VISP_HAVE_CATCH2)
42 #include <catch_amalgamated.hpp>
47 #include <visp3/core/vpGaussRand.h>
48 #include <visp3/core/vpHomogeneousMatrix.h>
49 #include <visp3/core/vpIoTools.h>
50 #include <visp3/core/vpMath.h>
51 #include <visp3/core/vpPoint.h>
52 #include <visp3/vision/vpPose.h>
54 #ifdef ENABLE_VISP_NAMESPACE
60 #if (VISP_HAVE_DATASET_VERSION >= 0x030300)
70 int checkInlierIndex(
const std::vector<unsigned int> &vectorOfFoundInlierIndex,
71 const std::vector<bool> &vectorOfOutlierFlags)
73 int nbInlierIndexOk = 0;
75 for (std::vector<unsigned int>::const_iterator it = vectorOfFoundInlierIndex.begin();
76 it != vectorOfFoundInlierIndex.end(); ++it) {
77 if (!vectorOfOutlierFlags[*it]) {
82 return nbInlierIndexOk;
85 bool checkInlierPoints(
const std::vector<vpPoint> &vectorOfFoundInlierPoints,
86 const std::vector<unsigned int> &vectorOfFoundInlierIndex,
87 const std::vector<vpPoint> &bunnyModelPoints_noisy)
89 for (
size_t i = 0; i < vectorOfFoundInlierPoints.size(); i++) {
90 if (!samePoints(vectorOfFoundInlierPoints[i], bunnyModelPoints_noisy[vectorOfFoundInlierIndex[i]])) {
91 std::cerr <<
"Problem with the inlier index and the corresponding "
94 std::cerr <<
"Returned inliers: oX=" << std::setprecision(std::numeric_limits<double>::max_digits10)
95 << vectorOfFoundInlierPoints[i].get_oX() <<
", oY=" << vectorOfFoundInlierPoints[i].get_oY()
96 <<
", oZ=" << vectorOfFoundInlierPoints[i].get_oZ() <<
" ; x=" << vectorOfFoundInlierPoints[i].get_x()
97 <<
", y=" << vectorOfFoundInlierPoints[i].get_y() << std::endl;
98 const vpPoint &pt = bunnyModelPoints_noisy[vectorOfFoundInlierIndex[i]];
99 std::cerr <<
"Object points: oX=" << std::setprecision(std::numeric_limits<double>::max_digits10) << pt.
get_oX()
109 void readBunnyModelPoints(
const std::string &filename, std::vector<vpPoint> &bunnyModelPoints,
110 std::vector<vpPoint> &bunnyModelPoints_noisy)
113 std::ifstream file(filename);
114 if (!file.is_open()) {
124 double oX = 0, oY = 0, oZ = 0;
126 while (file >> oX >> oY >> oZ) {
129 bunnyModelPoints.push_back(pt);
134 bunnyModelPoints_noisy.push_back(pt);
138 std::cout <<
"The raw model contains " << bunnyModelPoints.size() <<
" points." << std::endl;
139 std::cout <<
"cMo_groundTruth=\n" << cMo_groundTruth << std::endl << std::endl;
142 bool testRansac(
const std::vector<vpPoint> &bunnyModelPoints_original,
143 const std::vector<vpPoint> &bunnyModelPoints_noisy_original,
size_t nb_model_points,
144 bool test_duplicate,
bool test_degenerate)
146 std::vector<vpPoint> bunnyModelPoints = bunnyModelPoints_original;
147 std::vector<vpPoint> bunnyModelPoints_noisy = bunnyModelPoints_noisy_original;
149 if (nb_model_points > 0) {
150 bunnyModelPoints.resize(nb_model_points);
151 bunnyModelPoints_noisy.resize(nb_model_points);
154 vpPose ground_truth_pose, real_pose;
156 ground_truth_pose.
addPoints(bunnyModelPoints);
157 real_pose.
addPoints(bunnyModelPoints_noisy);
161 std::cout <<
"\ncMo estimated using VVS on data with small gaussian noise:\n" << cMo_estimated << std::endl;
162 std::cout <<
"Corresponding residual: " << r_vvs << std::endl;
164 size_t nbOutliers = (size_t)(0.35 * bunnyModelPoints_noisy.size());
167 std::vector<bool> vectorOfOutlierFlags(bunnyModelPoints_noisy.size(),
false);
169 for (
size_t i = 0; i < nbOutliers; i++) {
170 bunnyModelPoints_noisy[i].set_x(bunnyModelPoints_noisy[i].get_x() + noise());
171 bunnyModelPoints_noisy[i].set_y(bunnyModelPoints_noisy[i].get_y() + noise());
172 vectorOfOutlierFlags[i] =
true;
175 if (test_duplicate) {
177 size_t nbDuplicatePoints = 100;
178 for (
size_t i = 0; i < nbDuplicatePoints; i++) {
179 size_t index = (size_t)rand() % bunnyModelPoints_noisy.size();
180 vpPoint duplicatePoint = bunnyModelPoints_noisy[index];
181 bunnyModelPoints_noisy.push_back(duplicatePoint);
182 vectorOfOutlierFlags.push_back(
true);
186 if (test_degenerate) {
188 size_t nbDegeneratePoints = 100;
189 double degenerate_tolerence = 9.999e-7;
192 std::vector<vpPoint> listOfDegeneratePoints;
193 for (
size_t i = 0; i < nbDegeneratePoints; i++) {
194 size_t index = (size_t)rand() % bunnyModelPoints_noisy.size();
195 vpPoint degeneratePoint = bunnyModelPoints_noisy[index];
198 degeneratePoint.
set_oX(degeneratePoint.
get_oX() + degenerate_tolerence);
199 degeneratePoint.
set_oY(degeneratePoint.
get_oY() + degenerate_tolerence);
200 degeneratePoint.
set_oZ(degeneratePoint.
get_oZ() - degenerate_tolerence);
203 listOfDegeneratePoints.push_back(degeneratePoint);
206 index = (size_t)rand() % bunnyModelPoints_noisy.size();
207 degeneratePoint = bunnyModelPoints_noisy[index];
209 degeneratePoint.
set_x(degeneratePoint.
get_x() + degenerate_tolerence);
210 degeneratePoint.
set_y(degeneratePoint.
get_y() - degenerate_tolerence);
213 listOfDegeneratePoints.push_back(degeneratePoint);
216 for (std::vector<vpPoint>::const_iterator it_degenerate = listOfDegeneratePoints.begin();
217 it_degenerate != listOfDegeneratePoints.end(); ++it_degenerate) {
218 bunnyModelPoints_noisy.push_back(*it_degenerate);
219 vectorOfOutlierFlags.push_back(
true);
224 std::vector<size_t> vectorOfIndex(bunnyModelPoints_noisy.size());
225 for (
size_t i = 0; i < vectorOfIndex.size(); i++) {
226 vectorOfIndex[i] = i;
230 std::random_device rng;
231 std::mt19937 urng(rng());
232 std::shuffle(vectorOfIndex.begin(), vectorOfIndex.end(), urng);
234 std::vector<vpPoint> bunnyModelPoints_noisy_tmp = bunnyModelPoints_noisy;
235 bunnyModelPoints_noisy.clear();
236 std::vector<bool> vectorOfOutlierFlags_tmp = vectorOfOutlierFlags;
237 vectorOfOutlierFlags.clear();
238 for (std::vector<size_t>::const_iterator it = vectorOfIndex.begin(); it != vectorOfIndex.end(); ++it) {
239 bunnyModelPoints_noisy.push_back(bunnyModelPoints_noisy_tmp[*it]);
240 vectorOfOutlierFlags.push_back(vectorOfOutlierFlags_tmp[*it]);
245 vpPose pose_ransac, pose_ransac2;
247 vpPose pose_ransac_parallel, pose_ransac_parallel2;
255 for (std::vector<vpPoint>::const_iterator it = bunnyModelPoints_noisy.begin(); it != bunnyModelPoints_noisy.end();
260 pose_ransac.
addPoints(bunnyModelPoints_noisy);
261 pose_ransac2.
addPoints(bunnyModelPoints_noisy);
262 pose_ransac_parallel.
addPoints(bunnyModelPoints_noisy);
263 pose_ransac_parallel2.
addPoints(bunnyModelPoints_noisy);
266 std::cout <<
"\nNumber of model points in the noisy data vector: " << bunnyModelPoints_noisy.size() <<
" points."
270 unsigned int nbInlierToReachConsensus = (
unsigned int)(60.0 * (
double)(bunnyModelPoints_noisy.size()) / 100.0);
271 double threshold = 0.001;
290 std::cout <<
"Number of RANSAC iterations to ensure p=0.99 and epsilon=0.4: " << ransac_iterations << std::endl;
294 chrono_RANSAC.
start();
296 chrono_RANSAC.
stop();
298 std::cout <<
"\ncMo estimated with RANSAC (1000 iterations) on noisy data:\n" << cMo_estimated_RANSAC << std::endl;
299 std::cout <<
"Computation time: " << chrono_RANSAC.
getDurationMs() <<
" ms" << std::endl;
301 double r_RANSAC_estimated = ground_truth_pose.
computeResidual(cMo_estimated_RANSAC);
302 std::cout <<
"Corresponding residual (1000 iterations): " << r_RANSAC_estimated << std::endl;
305 chrono_RANSAC.
start();
307 chrono_RANSAC.
stop();
309 std::cout <<
"\ncMo estimated with RANSAC (" << ransac_iterations <<
" iterations) on noisy data:\n"
310 << cMo_estimated_RANSAC_2 << std::endl;
311 std::cout <<
"Computation time: " << chrono_RANSAC.
getDurationMs() <<
" ms" << std::endl;
313 double r_RANSAC_estimated_2 = ground_truth_pose.
computeResidual(cMo_estimated_RANSAC_2);
314 std::cout <<
"Corresponding residual (" << ransac_iterations <<
" iterations): " << r_RANSAC_estimated_2 << std::endl;
317 std::cout <<
"\ncMo estimated with only VVS on noisy data:\n" << cMo_estimated << std::endl;
320 std::cout <<
"Corresponding residual: " << r_estimated << std::endl;
324 chrono_RANSAC_parallel.
start();
326 chrono_RANSAC_parallel.
stop();
328 std::cout <<
"\ncMo estimated with parallel RANSAC (1000 iterations) on "
330 << cMo_estimated_RANSAC_parallel << std::endl;
331 std::cout <<
"Computation time: " << chrono_RANSAC_parallel.
getDurationMs() <<
" ms" << std::endl;
333 double r_RANSAC_estimated_parallel = ground_truth_pose.
computeResidual(cMo_estimated_RANSAC_parallel);
334 std::cout <<
"Corresponding residual (1000 iterations): " << r_RANSAC_estimated_parallel << std::endl;
338 chrono_RANSAC_parallel2.
start();
340 chrono_RANSAC_parallel2.
stop();
342 std::cout <<
"\ncMo estimated with parallel RANSAC (" << ransac_iterations <<
" iterations) on noisy data:\n"
343 << cMo_estimated_RANSAC_parallel2 << std::endl;
344 std::cout <<
"Computation time: " << chrono_RANSAC_parallel2.
getDurationMs() <<
" ms" << std::endl;
346 double r_RANSAC_estimated_parallel2 = ground_truth_pose.
computeResidual(cMo_estimated_RANSAC_parallel2);
347 std::cout <<
"Corresponding residual (" << ransac_iterations <<
" iterations): " << r_RANSAC_estimated_parallel2
352 int nbInlierIndexOk = checkInlierIndex(vectorOfFoundInlierIndex, vectorOfOutlierFlags);
354 int nbTrueInlierIndex = (int)std::count(vectorOfOutlierFlags.begin(), vectorOfOutlierFlags.end(),
false);
355 std::cout <<
"\nThere are " << nbInlierIndexOk <<
" true inliers found, " << vectorOfFoundInlierIndex.size()
356 <<
" inliers returned and " << nbTrueInlierIndex <<
" true inliers." << std::endl;
359 std::vector<vpPoint> vectorOfFoundInlierPoints = pose_ransac.
getRansacInliers();
361 if (vectorOfFoundInlierPoints.size() != vectorOfFoundInlierIndex.size()) {
362 std::cerr <<
"The number of inlier index is different from the number of "
367 if (!checkInlierPoints(vectorOfFoundInlierPoints, vectorOfFoundInlierIndex, bunnyModelPoints_noisy)) {
373 std::cout <<
"\nCheck for RANSAC iterations: " << ransac_iterations << std::endl;
375 nbInlierIndexOk = checkInlierIndex(vectorOfFoundInlierIndex_2, vectorOfOutlierFlags);
377 std::cout <<
"There are " << nbInlierIndexOk <<
" true inliers found, " << vectorOfFoundInlierIndex_2.size()
378 <<
" inliers returned and " << nbTrueInlierIndex <<
" true inliers." << std::endl;
381 std::vector<vpPoint> vectorOfFoundInlierPoints_2 = pose_ransac2.
getRansacInliers();
382 if (vectorOfFoundInlierPoints_2.size() != vectorOfFoundInlierIndex_2.size()) {
383 std::cerr <<
"The number of inlier index is different from the number of "
388 if (!checkInlierPoints(vectorOfFoundInlierPoints_2, vectorOfFoundInlierIndex_2, bunnyModelPoints_noisy)) {
394 std::cout <<
"\nCheck for parallel RANSAC (1000 iterations)" << std::endl;
395 std::vector<unsigned int> vectorOfFoundInlierIndex_parallel = pose_ransac_parallel.
getRansacInlierIndex();
396 nbInlierIndexOk = checkInlierIndex(vectorOfFoundInlierIndex_parallel, vectorOfOutlierFlags);
398 std::cout <<
"There are " << nbInlierIndexOk <<
" true inliers found, " << vectorOfFoundInlierIndex_parallel.size()
399 <<
" inliers returned and " << nbTrueInlierIndex <<
" true inliers." << std::endl;
402 std::vector<vpPoint> vectorOfFoundInlierPoints_parallel = pose_ransac_parallel.
getRansacInliers();
403 if (vectorOfFoundInlierPoints_parallel.size() != vectorOfFoundInlierIndex_parallel.size()) {
404 std::cerr <<
"The number of inlier index is different from the number "
409 if (!checkInlierPoints(vectorOfFoundInlierPoints_parallel, vectorOfFoundInlierIndex_parallel,
410 bunnyModelPoints_noisy)) {
416 std::cout <<
"\nCheck for parallel RANSAC (" << ransac_iterations <<
" iterations)" << std::endl;
417 std::vector<unsigned int> vectorOfFoundInlierIndex_parallel2 = pose_ransac_parallel2.
getRansacInlierIndex();
418 nbInlierIndexOk = checkInlierIndex(vectorOfFoundInlierIndex_parallel2, vectorOfOutlierFlags);
420 std::cout <<
"There are " << nbInlierIndexOk <<
" true inliers found, " << vectorOfFoundInlierIndex_parallel2.size()
421 <<
" inliers returned and " << nbTrueInlierIndex <<
" true inliers." << std::endl;
424 std::vector<vpPoint> vectorOfFoundInlierPoints_parallel2 = pose_ransac_parallel2.
getRansacInliers();
425 if (vectorOfFoundInlierPoints_parallel2.size() != vectorOfFoundInlierIndex_parallel2.size()) {
426 std::cerr <<
"The number of inlier index is different from the number "
431 if (!checkInlierPoints(vectorOfFoundInlierPoints_parallel2, vectorOfFoundInlierIndex_parallel2,
432 bunnyModelPoints_noisy)) {
436 if (r_RANSAC_estimated > threshold ) {
437 std::cerr <<
"The pose estimated with the RANSAC method is badly estimated!" << std::endl;
438 std::cerr <<
"r_RANSAC_estimated=" << r_RANSAC_estimated << std::endl;
439 std::cerr <<
"threshold=" << threshold << std::endl;
443 if (r_RANSAC_estimated_parallel > threshold) {
444 std::cerr <<
"The pose estimated with the parallel RANSAC method is "
447 std::cerr <<
"r_RANSAC_estimated_parallel=" << r_RANSAC_estimated_parallel << std::endl;
448 std::cerr <<
"threshold=" << threshold << std::endl;
451 std::cout <<
"The pose estimated with the RANSAC method is well estimated!" << std::endl;
459 TEST_CASE(
"Print RANSAC number of iterations",
"[ransac_pose]")
461 const int sample_sizes[] = { 2, 3, 4, 5, 6, 7, 8 };
462 const double epsilon[] = { 0.05, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5 };
465 const std::string spacing =
" ";
467 std::cout << spacing <<
" outliers percentage\n"
469 for (
int cpt2 = 0; cpt2 < 7; cpt2++) {
470 std::cout << std::setfill(
' ') << std::setw(5) << epsilon[cpt2] <<
" ";
472 std::cout << std::endl;
474 std::cout << std::setfill(
' ') << std::setw(7) <<
"+";
475 for (
int cpt2 = 0; cpt2 < 6; cpt2++) {
476 std::cout << std::setw(7) <<
"-------";
478 std::cout << std::endl;
480 for (
int cpt1 = 0; cpt1 < 7; cpt1++) {
481 std::cout << std::setfill(
' ') << std::setw(6) << sample_sizes[cpt1] <<
"|";
483 for (
int cpt2 = 0; cpt2 < 7; cpt2++) {
485 std::cout << std::setfill(
' ') << std::setw(6) << ransac_iters;
487 std::cout << std::endl;
489 std::cout << std::endl;
492 #if (VISP_HAVE_DATASET_VERSION >= 0x030300)
493 TEST_CASE(
"RANSAC pose estimation tests",
"[ransac_pose]")
495 const std::vector<size_t> model_sizes = { 10, 20, 50, 100, 200, 500, 1000, 0, 0 };
496 const std::vector<bool> duplicates = {
false,
false,
false,
false,
false,
false,
false,
false,
true };
497 const std::vector<bool> degenerates = {
false,
false,
false,
false,
false,
false,
true,
true,
true };
503 std::vector<vpPoint> bunnyModelPoints, bunnyModelPoints_noisy_original;
504 readBunnyModelPoints(model_filename, bunnyModelPoints, bunnyModelPoints_noisy_original);
505 CHECK(bunnyModelPoints.size() == bunnyModelPoints_noisy_original.size());
507 for (
size_t i = 0; i < model_sizes.size(); i++) {
508 std::cout <<
"\n\n===============================================================================" << std::endl;
509 if (model_sizes[i] == 0) {
510 std::cout <<
"Test on " << bunnyModelPoints_noisy_original.size() <<
" model points." << std::endl;
513 std::cout <<
"Test on " << model_sizes[i] <<
" model points." << std::endl;
515 std::cout <<
"Test duplicate: " << duplicates[i] <<
" ; Test degenerate: " << degenerates[i] << std::endl;
517 CHECK(testRansac(bunnyModelPoints, bunnyModelPoints_noisy_original, model_sizes[i], duplicates[i], degenerates[i]));
522 int main(
int argc,
char *argv[])
524 #if defined(__mips__) || defined(__mips) || defined(mips) || defined(__MIPS__)
529 #if (defined(VISP_HAVE_LAPACK) || defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_OPENCV))
531 Catch::Session session;
534 session.applyCommandLine(argc, argv);
536 int numFailed = session.run();
543 std::cout <<
"Cannot run this example: install Lapack, Eigen3 or OpenCV" << std::endl;
548 int main() {
return EXIT_SUCCESS; }
void start(bool reset=true)
Class for generating random number with normal probability density.
Implementation of an homogeneous matrix and operations on such kind of matrices.
static double rad(double deg)
static bool equal(double x, double y, double threshold=0.001)
Class that defines a 3D point in the object frame and allows forward projection of a 3D point in the ...
double get_oX() const
Get the point oX coordinate in the object frame.
void set_x(double x)
Set the point x coordinate in the image plane.
double get_y() const
Get the point y coordinate in the image plane.
double get_oZ() const
Get the point oZ coordinate in the object frame.
void set_oY(double oY)
Set the point oY coordinate in the object frame.
double get_x() const
Get the point x coordinate in the image plane.
void set_oZ(double oZ)
Set the point oZ coordinate in the object frame.
void set_oX(double oX)
Set the point oX coordinate in the object frame.
double get_oY() const
Get the point oY coordinate in the object frame.
void set_y(double y)
Set the point y coordinate in the image plane.
Class used for pose computation from N points (pose from point only). Some of the algorithms implemen...
void setRansacMaxTrials(const int &rM)
static int computeRansacIterations(double probability, double epsilon, const int sampleSize=4, int maxIterations=2000)
void addPoint(const vpPoint &P)
void setRansacNbInliersToReachConsensus(const unsigned int &nbC)
@ DEMENTHON_LAGRANGE_VIRTUAL_VS
bool computePose(vpPoseMethodType method, vpHomogeneousMatrix &cMo, FuncCheckValidityPose func=nullptr)
std::vector< unsigned int > getRansacInlierIndex() const
void addPoints(const std::vector< vpPoint > &lP)
double computeResidual(const vpHomogeneousMatrix &cMo) const
Compute and return the sum of squared residuals expressed in meter^2 for the pose matrix cMo.
void setRansacFilterFlag(const RANSAC_FILTER_FLAGS &flag)
@ PREFILTER_DEGENERATE_POINTS
void setUseParallelRansac(bool use)
std::vector< vpPoint > getRansacInliers() const
void setRansacThreshold(const double &t)
Implementation of a rotation vector as Euler angle minimal representation.
Implementation of a rotation vector as axis-angle minimal representation.
Class that consider the case of a translation vector.