#include <visp3/core/vpConfig.h>
#if defined(VISP_HAVE_CATCH2)
#include <catch_amalgamated.hpp>
#include <algorithm>
#include <iomanip>
#include <map>
#include <visp3/core/vpGaussRand.h>
#include <visp3/core/vpHomogeneousMatrix.h>
#include <visp3/core/vpIoTools.h>
#include <visp3/core/vpMath.h>
#include <visp3/core/vpPoint.h>
#include <visp3/vision/vpPose.h>
#ifdef ENABLE_VISP_NAMESPACE
#endif
namespace
{
#if (VISP_HAVE_DATASET_VERSION >= 0x030300)
{
}
int checkInlierIndex(const std::vector<unsigned int> &vectorOfFoundInlierIndex,
const std::vector<bool> &vectorOfOutlierFlags)
{
int nbInlierIndexOk = 0;
for (std::vector<unsigned int>::const_iterator it = vectorOfFoundInlierIndex.begin();
it != vectorOfFoundInlierIndex.end(); ++it) {
if (!vectorOfOutlierFlags[*it]) {
nbInlierIndexOk++;
}
}
return nbInlierIndexOk;
}
bool checkInlierPoints(const std::vector<vpPoint> &vectorOfFoundInlierPoints,
const std::vector<unsigned int> &vectorOfFoundInlierIndex,
const std::vector<vpPoint> &bunnyModelPoints_noisy)
{
for (size_t i = 0; i < vectorOfFoundInlierPoints.size(); i++) {
if (!samePoints(vectorOfFoundInlierPoints[i], bunnyModelPoints_noisy[vectorOfFoundInlierIndex[i]])) {
std::cerr << "Problem with the inlier index and the corresponding "
"inlier point!"
<< std::endl;
std::cerr << "Returned inliers: oX=" << std::setprecision(std::numeric_limits<double>::max_digits10)
<< vectorOfFoundInlierPoints[i].get_oX() << ", oY=" << vectorOfFoundInlierPoints[i].get_oY()
<< ", oZ=" << vectorOfFoundInlierPoints[i].get_oZ() << " ; x=" << vectorOfFoundInlierPoints[i].get_x()
<< ", y=" << vectorOfFoundInlierPoints[i].get_y() << std::endl;
const vpPoint &pt = bunnyModelPoints_noisy[vectorOfFoundInlierIndex[i]];
std::cerr <<
"Object points: oX=" << std::setprecision(std::numeric_limits<double>::max_digits10) << pt.
get_oX()
<< std::endl;
return false;
}
}
return true;
}
void readBunnyModelPoints(const std::string &filename, std::vector<vpPoint> &bunnyModelPoints,
std::vector<vpPoint> &bunnyModelPoints_noisy)
{
std::ifstream file(filename);
if (!file.is_open()) {
return;
}
double oX = 0, oY = 0, oZ = 0;
while (file >> oX >> oY >> oZ) {
bunnyModelPoints.push_back(pt);
bunnyModelPoints_noisy.push_back(pt);
}
std::cout << "The raw model contains " << bunnyModelPoints.size() << " points." << std::endl;
std::cout << "cMo_groundTruth=\n" << cMo_groundTruth << std::endl << std::endl;
}
bool testRansac(const std::vector<vpPoint> &bunnyModelPoints_original,
const std::vector<vpPoint> &bunnyModelPoints_noisy_original, size_t nb_model_points,
bool test_duplicate, bool test_degenerate)
{
std::vector<vpPoint> bunnyModelPoints = bunnyModelPoints_original;
std::vector<vpPoint> bunnyModelPoints_noisy = bunnyModelPoints_noisy_original;
if (nb_model_points > 0) {
bunnyModelPoints.resize(nb_model_points);
bunnyModelPoints_noisy.resize(nb_model_points);
}
vpPose ground_truth_pose, real_pose;
ground_truth_pose.
addPoints(bunnyModelPoints);
std::cout << "\ncMo estimated using VVS on data with small gaussian noise:\n" << cMo_estimated << std::endl;
std::cout << "Corresponding residual: " << r_vvs << std::endl;
size_t nbOutliers = static_cast<size_t>(0.35 * bunnyModelPoints_noisy.size());
std::vector<bool> vectorOfOutlierFlags(bunnyModelPoints_noisy.size(), false);
for (size_t i = 0; i < nbOutliers; i++) {
bunnyModelPoints_noisy[i].set_x(bunnyModelPoints_noisy[i].get_x() + noise());
bunnyModelPoints_noisy[i].set_y(bunnyModelPoints_noisy[i].get_y() + noise());
vectorOfOutlierFlags[i] = true;
}
if (test_duplicate) {
size_t nbDuplicatePoints = 100;
for (size_t i = 0; i < nbDuplicatePoints; i++) {
size_t index = static_cast<size_t>(rand()) % bunnyModelPoints_noisy.size();
vpPoint duplicatePoint = bunnyModelPoints_noisy[index];
bunnyModelPoints_noisy.push_back(duplicatePoint);
vectorOfOutlierFlags.push_back(true);
}
}
if (test_degenerate) {
size_t nbDegeneratePoints = 100;
double degenerate_tolerence = 9.999e-7;
std::vector<vpPoint> listOfDegeneratePoints;
for (size_t i = 0; i < nbDegeneratePoints; i++) {
size_t index = static_cast<size_t>(rand()) % bunnyModelPoints_noisy.size();
vpPoint degeneratePoint = bunnyModelPoints_noisy[index];
degeneratePoint.
set_oX(degeneratePoint.
get_oX() + degenerate_tolerence);
degeneratePoint.
set_oY(degeneratePoint.
get_oY() + degenerate_tolerence);
degeneratePoint.
set_oZ(degeneratePoint.
get_oZ() - degenerate_tolerence);
listOfDegeneratePoints.push_back(degeneratePoint);
index = static_cast<size_t>(rand()) % bunnyModelPoints_noisy.size();
degeneratePoint = bunnyModelPoints_noisy[index];
degeneratePoint.
set_x(degeneratePoint.
get_x() + degenerate_tolerence);
degeneratePoint.
set_y(degeneratePoint.
get_y() - degenerate_tolerence);
listOfDegeneratePoints.push_back(degeneratePoint);
}
for (std::vector<vpPoint>::const_iterator it_degenerate = listOfDegeneratePoints.begin();
it_degenerate != listOfDegeneratePoints.end(); ++it_degenerate) {
bunnyModelPoints_noisy.push_back(*it_degenerate);
vectorOfOutlierFlags.push_back(true);
}
}
std::vector<size_t> vectorOfIndex(bunnyModelPoints_noisy.size());
for (size_t i = 0; i < vectorOfIndex.size(); i++) {
vectorOfIndex[i] = i;
}
std::random_device rng;
std::mt19937 urng(rng());
std::shuffle(vectorOfIndex.begin(), vectorOfIndex.end(), urng);
std::vector<vpPoint> bunnyModelPoints_noisy_tmp = bunnyModelPoints_noisy;
bunnyModelPoints_noisy.clear();
std::vector<bool> vectorOfOutlierFlags_tmp = vectorOfOutlierFlags;
vectorOfOutlierFlags.clear();
for (std::vector<size_t>::const_iterator it = vectorOfIndex.begin(); it != vectorOfIndex.end(); ++it) {
bunnyModelPoints_noisy.push_back(bunnyModelPoints_noisy_tmp[*it]);
vectorOfOutlierFlags.push_back(vectorOfOutlierFlags_tmp[*it]);
}
vpPose pose_ransac, pose_ransac2;
vpPose pose_ransac_parallel, pose_ransac_parallel2;
for (std::vector<vpPoint>::const_iterator it = bunnyModelPoints_noisy.begin(); it != bunnyModelPoints_noisy.end();
++it) {
}
pose_ransac.
addPoints(bunnyModelPoints_noisy);
pose_ransac2.
addPoints(bunnyModelPoints_noisy);
pose_ransac_parallel.
addPoints(bunnyModelPoints_noisy);
pose_ransac_parallel2.
addPoints(bunnyModelPoints_noisy);
std::cout << "\nNumber of model points in the noisy data vector: " << bunnyModelPoints_noisy.size() << " points."
<< std::endl
<< std::endl;
unsigned int nbInlierToReachConsensus = static_cast<unsigned int>(60.0 * static_cast<double>(bunnyModelPoints_noisy.size()) / 100.0);
double threshold = 0.001;
std::cout << "Number of RANSAC iterations to ensure p=0.99 and epsilon=0.4: " << ransac_iterations << std::endl;
std::cout << "\ncMo estimated with RANSAC (1000 iterations) on noisy data:\n" << cMo_estimated_RANSAC << std::endl;
std::cout <<
"Computation time: " << chrono_RANSAC.
getDurationMs() <<
" ms" << std::endl;
double r_RANSAC_estimated = ground_truth_pose.
computeResidual(cMo_estimated_RANSAC);
std::cout << "Corresponding residual (1000 iterations): " << r_RANSAC_estimated << std::endl;
std::cout << "\ncMo estimated with RANSAC (" << ransac_iterations << " iterations) on noisy data:\n"
<< cMo_estimated_RANSAC_2 << std::endl;
std::cout <<
"Computation time: " << chrono_RANSAC.
getDurationMs() <<
" ms" << std::endl;
double r_RANSAC_estimated_2 = ground_truth_pose.
computeResidual(cMo_estimated_RANSAC_2);
std::cout << "Corresponding residual (" << ransac_iterations << " iterations): " << r_RANSAC_estimated_2 << std::endl;
std::cout << "\ncMo estimated with only VVS on noisy data:\n" << cMo_estimated << std::endl;
std::cout << "Corresponding residual: " << r_estimated << std::endl;
chrono_RANSAC_parallel.
start();
chrono_RANSAC_parallel.
stop();
std::cout << "\ncMo estimated with parallel RANSAC (1000 iterations) on "
"noisy data:\n"
<< cMo_estimated_RANSAC_parallel << std::endl;
std::cout <<
"Computation time: " << chrono_RANSAC_parallel.
getDurationMs() <<
" ms" << std::endl;
double r_RANSAC_estimated_parallel = ground_truth_pose.
computeResidual(cMo_estimated_RANSAC_parallel);
std::cout << "Corresponding residual (1000 iterations): " << r_RANSAC_estimated_parallel << std::endl;
chrono_RANSAC_parallel2.
start();
chrono_RANSAC_parallel2.
stop();
std::cout << "\ncMo estimated with parallel RANSAC (" << ransac_iterations << " iterations) on noisy data:\n"
<< cMo_estimated_RANSAC_parallel2 << std::endl;
std::cout <<
"Computation time: " << chrono_RANSAC_parallel2.
getDurationMs() <<
" ms" << std::endl;
double r_RANSAC_estimated_parallel2 = ground_truth_pose.
computeResidual(cMo_estimated_RANSAC_parallel2);
std::cout << "Corresponding residual (" << ransac_iterations << " iterations): " << r_RANSAC_estimated_parallel2
<< std::endl;
int nbInlierIndexOk = checkInlierIndex(vectorOfFoundInlierIndex, vectorOfOutlierFlags);
int nbTrueInlierIndex = static_cast<int>(std::count(vectorOfOutlierFlags.begin(), vectorOfOutlierFlags.end(), false));
std::cout << "\nThere are " << nbInlierIndexOk << " true inliers found, " << vectorOfFoundInlierIndex.size()
<< " inliers returned and " << nbTrueInlierIndex << " true inliers." << std::endl;
std::vector<vpPoint> vectorOfFoundInlierPoints = pose_ransac.
getRansacInliers();
if (vectorOfFoundInlierPoints.size() != vectorOfFoundInlierIndex.size()) {
std::cerr << "The number of inlier index is different from the number of "
"inlier points!"
<< std::endl;
return false;
}
if (!checkInlierPoints(vectorOfFoundInlierPoints, vectorOfFoundInlierIndex, bunnyModelPoints_noisy)) {
return false;
}
std::cout << "\nCheck for RANSAC iterations: " << ransac_iterations << std::endl;
nbInlierIndexOk = checkInlierIndex(vectorOfFoundInlierIndex_2, vectorOfOutlierFlags);
std::cout << "There are " << nbInlierIndexOk << " true inliers found, " << vectorOfFoundInlierIndex_2.size()
<< " inliers returned and " << nbTrueInlierIndex << " true inliers." << std::endl;
std::vector<vpPoint> vectorOfFoundInlierPoints_2 = pose_ransac2.
getRansacInliers();
if (vectorOfFoundInlierPoints_2.size() != vectorOfFoundInlierIndex_2.size()) {
std::cerr << "The number of inlier index is different from the number of "
"inlier points!"
<< std::endl;
return false;
}
if (!checkInlierPoints(vectorOfFoundInlierPoints_2, vectorOfFoundInlierIndex_2, bunnyModelPoints_noisy)) {
return false;
}
std::cout << "\nCheck for parallel RANSAC (1000 iterations)" << std::endl;
std::vector<unsigned int> vectorOfFoundInlierIndex_parallel = pose_ransac_parallel.
getRansacInlierIndex();
nbInlierIndexOk = checkInlierIndex(vectorOfFoundInlierIndex_parallel, vectorOfOutlierFlags);
std::cout << "There are " << nbInlierIndexOk << " true inliers found, " << vectorOfFoundInlierIndex_parallel.size()
<< " inliers returned and " << nbTrueInlierIndex << " true inliers." << std::endl;
std::vector<vpPoint> vectorOfFoundInlierPoints_parallel = pose_ransac_parallel.
getRansacInliers();
if (vectorOfFoundInlierPoints_parallel.size() != vectorOfFoundInlierIndex_parallel.size()) {
std::cerr << "The number of inlier index is different from the number "
"of inlier points!"
<< std::endl;
return false;
}
if (!checkInlierPoints(vectorOfFoundInlierPoints_parallel, vectorOfFoundInlierIndex_parallel,
bunnyModelPoints_noisy)) {
return false;
}
std::cout << "\nCheck for parallel RANSAC (" << ransac_iterations << " iterations)" << std::endl;
std::vector<unsigned int> vectorOfFoundInlierIndex_parallel2 = pose_ransac_parallel2.
getRansacInlierIndex();
nbInlierIndexOk = checkInlierIndex(vectorOfFoundInlierIndex_parallel2, vectorOfOutlierFlags);
std::cout << "There are " << nbInlierIndexOk << " true inliers found, " << vectorOfFoundInlierIndex_parallel2.size()
<< " inliers returned and " << nbTrueInlierIndex << " true inliers." << std::endl;
std::vector<vpPoint> vectorOfFoundInlierPoints_parallel2 = pose_ransac_parallel2.
getRansacInliers();
if (vectorOfFoundInlierPoints_parallel2.size() != vectorOfFoundInlierIndex_parallel2.size()) {
std::cerr << "The number of inlier index is different from the number "
"of inlier points!"
<< std::endl;
return false;
}
if (!checkInlierPoints(vectorOfFoundInlierPoints_parallel2, vectorOfFoundInlierIndex_parallel2,
bunnyModelPoints_noisy)) {
return false;
}
if (r_RANSAC_estimated > threshold ) {
std::cerr << "The pose estimated with the RANSAC method is badly estimated!" << std::endl;
std::cerr << "r_RANSAC_estimated=" << r_RANSAC_estimated << std::endl;
std::cerr << "threshold=" << threshold << std::endl;
return false;
}
else {
if (r_RANSAC_estimated_parallel > threshold) {
std::cerr << "The pose estimated with the parallel RANSAC method is "
"badly estimated!"
<< std::endl;
std::cerr << "r_RANSAC_estimated_parallel=" << r_RANSAC_estimated_parallel << std::endl;
std::cerr << "threshold=" << threshold << std::endl;
return false;
}
std::cout << "The pose estimated with the RANSAC method is well estimated!" << std::endl;
}
return true;
}
#endif
}
TEST_CASE("Print RANSAC number of iterations", "[ransac_pose]")
{
const int sample_sizes[] = { 2, 3, 4, 5, 6, 7, 8 };
const double epsilon[] = { 0.05, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5 };
const std::string spacing = " ";
std::cout << spacing << " outliers percentage\n"
<< "nb pts\\";
for (int cpt2 = 0; cpt2 < 7; cpt2++) {
std::cout << std::setfill(' ') << std::setw(5) << epsilon[cpt2] << " ";
}
std::cout << std::endl;
std::cout << std::setfill(' ') << std::setw(7) << "+";
for (int cpt2 = 0; cpt2 < 6; cpt2++) {
std::cout << std::setw(7) << "-------";
}
std::cout << std::endl;
for (int cpt1 = 0; cpt1 < 7; cpt1++) {
std::cout << std::setfill(' ') << std::setw(6) << sample_sizes[cpt1] << "|";
for (int cpt2 = 0; cpt2 < 7; cpt2++) {
std::cout << std::setfill(' ') << std::setw(6) << ransac_iters;
}
std::cout << std::endl;
}
std::cout << std::endl;
}
#if (VISP_HAVE_DATASET_VERSION >= 0x030300)
TEST_CASE("RANSAC pose estimation tests", "[ransac_pose]")
{
const std::vector<size_t> model_sizes = { 10, 20, 50, 100, 200, 500, 1000, 0, 0 };
const std::vector<bool> duplicates = { false, false, false, false, false, false, false, false, true };
const std::vector<bool> degenerates = { false, false, false, false, false, false, true, true, true };
std::vector<vpPoint> bunnyModelPoints, bunnyModelPoints_noisy_original;
readBunnyModelPoints(model_filename, bunnyModelPoints, bunnyModelPoints_noisy_original);
CHECK(bunnyModelPoints.size() == bunnyModelPoints_noisy_original.size());
for (
size_t i = 0;
i < model_sizes.size();
i++) {
std::cout << "\n\n===============================================================================" << std::endl;
if (model_sizes[i] == 0) {
std::cout << "Test on " << bunnyModelPoints_noisy_original.size() << " model points." << std::endl;
}
else {
std::cout <<
"Test on " << model_sizes[
i] <<
" model points." << std::endl;
}
std::cout <<
"Test duplicate: " << duplicates[
i] <<
" ; Test degenerate: " << degenerates[
i] << std::endl;
CHECK(testRansac(bunnyModelPoints, bunnyModelPoints_noisy_original, model_sizes[i], duplicates[i], degenerates[i]));
}
}
#endif
int main(int argc, char *argv[])
{
#if defined(__mips__) || defined(__mips) || defined(mips) || defined(__MIPS__)
return EXIT_SUCCESS;
#endif
#if (defined(VISP_HAVE_LAPACK) || defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_OPENCV))
Catch::Session session;
session.applyCommandLine(argc, argv);
int numFailed = session.run();
return numFailed;
#else
std::cout << "Cannot run this example: install Lapack, Eigen3 or OpenCV" << std::endl;
return EXIT_SUCCESS;
#endif
}
#else
int main() { return EXIT_SUCCESS; }
#endif
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)
std::vector< unsigned int > getRansacInlierIndex() const
@ DEMENTHON_LAGRANGE_VIRTUAL_VS
bool computePose(vpPoseMethodType method, vpHomogeneousMatrix &cMo, FuncCheckValidityPose func=nullptr)
std::vector< vpPoint > getRansacInliers() 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)
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.