39 #include <pcl/conversions.h> 40 #include <pcl/common/pca.h> 41 #include <pcl/common/io.h> 43 #ifdef RTABMAP_POINTMATCHER 45 #include "pointmatcher/PointMatcher.h" 46 #include "nabo/nabo.h" 47 typedef PointMatcher<float> PM;
48 typedef PM::DataPoints DP;
50 DP pclToDP(
const pcl::PointCloud<pcl::PointXYZ>::Ptr & pclCloud,
bool is2D)
53 typedef DP::Label Label;
54 typedef DP::Labels Labels;
55 typedef DP::View View;
57 if (pclCloud->empty())
65 std::vector<bool> isFeature;
66 featLabels.push_back(Label(
"x", 1));
67 isFeature.push_back(
true);
68 featLabels.push_back(Label(
"y", 1));
69 isFeature.push_back(
true);
72 featLabels.push_back(Label(
"z", 1));
73 isFeature.push_back(
true);
75 featLabels.push_back(Label(
"pad", 1));
78 DP cloud(featLabels, descLabels, pclCloud->size());
79 cloud.getFeatureViewByName(
"pad").setConstant(1);
82 View view(cloud.getFeatureViewByName(
"x"));
83 for(
unsigned int i=0; i<pclCloud->size(); ++i)
85 view(0, i) = pclCloud->at(i).x;
86 view(1, i) = pclCloud->at(i).y;
89 view(2, i) = pclCloud->at(i).z;
96 DP pclToDP(
const pcl::PointCloud<pcl::PointNormal>::Ptr & pclCloud,
bool is2D)
99 typedef DP::Label Label;
100 typedef DP::Labels Labels;
101 typedef DP::View View;
103 if (pclCloud->empty())
111 std::vector<bool> isFeature;
112 featLabels.push_back(Label(
"x", 1));
113 isFeature.push_back(
true);
114 featLabels.push_back(Label(
"y", 1));
115 isFeature.push_back(
true);
118 featLabels.push_back(Label(
"z", 1));
119 isFeature.push_back(
true);
122 descLabels.push_back(Label(
"normals", 3));
123 isFeature.push_back(
false);
124 isFeature.push_back(
false);
125 isFeature.push_back(
false);
127 featLabels.push_back(Label(
"pad", 1));
130 DP cloud(featLabels, descLabels, pclCloud->size());
131 cloud.getFeatureViewByName(
"pad").setConstant(1);
134 View view(cloud.getFeatureViewByName(
"x"));
135 View viewNormalX(cloud.getDescriptorRowViewByName(
"normals",0));
136 View viewNormalY(cloud.getDescriptorRowViewByName(
"normals",1));
137 View viewNormalZ(cloud.getDescriptorRowViewByName(
"normals",2));
138 for(
unsigned int i=0; i<pclCloud->size(); ++i)
140 view(0, i) = pclCloud->at(i).x;
141 view(1, i) = pclCloud->at(i).y;
144 view(2, i) = pclCloud->at(i).z;
146 viewNormalX(0, i) = pclCloud->at(i).normal_x;
147 viewNormalY(0, i) = pclCloud->at(i).normal_y;
148 viewNormalZ(0, i) = pclCloud->at(i).normal_z;
157 typedef DP::Label Label;
158 typedef DP::Labels Labels;
159 typedef DP::View View;
169 featLabels.push_back(Label(
"x", 1));
170 featLabels.push_back(Label(
"y", 1));
173 featLabels.push_back(Label(
"z", 1));
175 featLabels.push_back(Label(
"pad", 1));
179 descLabels.push_back(Label(
"normals", 3));
183 descLabels.push_back(Label(
"intensity", 1));
188 DP cloud(featLabels, descLabels, scan.
size());
189 cloud.getFeatureViewByName(
"pad").setConstant(1);
197 View view(cloud.getFeatureViewByName(
"x"));
198 View viewNormalX(nx!=-1?cloud.getDescriptorRowViewByName(
"normals",0):view);
199 View viewNormalY(nx!=-1?cloud.getDescriptorRowViewByName(
"normals",1):view);
200 View viewNormalZ(nx!=-1?cloud.getDescriptorRowViewByName(
"normals",2):view);
201 View viewIntensity(offsetI!=-1?cloud.getDescriptorRowViewByName(
"intensity",0):view);
203 for(
int i=0; i<scan.
size(); ++i)
205 const float * ptr = scan.
data().ptr<
float>(0, i);
209 if(hasLocalTransform)
213 cv::Point3f pt(ptr[0], ptr[1], scan.
is2d()?0:ptr[2]);
223 viewIntensity(0, oi) = ptr[offsetI];
232 pt.z=scan.
is2d()?0:ptr[2];
243 viewNormalX(0, oi) = pt.normal_x;
244 viewNormalY(0, oi) = pt.normal_y;
245 viewNormalZ(0, oi) = pt.normal_z;
249 viewIntensity(0, oi) = ptr[offsetI];
256 UWARN(
"Ignoring point %d with invalid data: pos=%f %f %f, normal=%f %f %f", i, ptr[0], ptr[1], scan.
is2d()?0:ptr[3], ptr[nx], ptr[ny], ptr[nz]);
261 view(0, oi) = ptr[0];
262 view(1, oi) = ptr[1];
265 view(2, oi) = ptr[2];
269 viewNormalX(0, oi) = ptr[nx];
270 viewNormalY(0, oi) = ptr[ny];
271 viewNormalZ(0, oi) = ptr[nz];
275 viewIntensity(0, oi) = ptr[offsetI];
281 UWARN(
"Ignoring point %d with invalid data: pos=%f %f %f, normal=%f %f %f", i, ptr[0], ptr[1], scan.
is2d()?0:ptr[3], ptr[nx], ptr[ny], ptr[nz]);
286 UWARN(
"Ignoring point %d with invalid data: pos=%f %f %f", i, ptr[0], ptr[1], scan.
is2d()?0:ptr[3]);
290 if(oi != scan.
size())
292 cloud.conservativeResize(oi);
298 void pclFromDP(
const DP & cloud, pcl::PointCloud<pcl::PointXYZ> & pclCloud)
301 typedef DP::ConstView ConstView;
303 if (cloud.features.cols() == 0)
306 pclCloud.resize(cloud.features.cols());
307 pclCloud.is_dense =
true;
310 ConstView view(cloud.getFeatureViewByName(
"x"));
311 bool is3D = cloud.featureExists(
"z");
312 for(
unsigned int i=0; i<pclCloud.size(); ++i)
314 pclCloud.at(i).x = view(0, i);
315 pclCloud.at(i).y = view(1, i);
316 pclCloud.at(i).z = is3D?view(2, i):0;
320 void pclFromDP(
const DP & cloud, pcl::PointCloud<pcl::PointNormal> & pclCloud)
323 typedef DP::ConstView ConstView;
325 if (cloud.features.cols() == 0)
328 pclCloud.resize(cloud.features.cols());
329 pclCloud.is_dense =
true;
332 ConstView view(cloud.getFeatureViewByName(
"x"));
333 bool is3D = cloud.featureExists(
"z");
334 ConstView viewNormalX(cloud.getDescriptorRowViewByName(
"normals",0));
335 ConstView viewNormalY(cloud.getDescriptorRowViewByName(
"normals",1));
336 ConstView viewNormalZ(cloud.getDescriptorRowViewByName(
"normals",2));
337 for(
unsigned int i=0; i<pclCloud.size(); ++i)
339 pclCloud.at(i).x = view(0, i);
340 pclCloud.at(i).y = view(1, i);
341 pclCloud.at(i).z = is3D?view(2, i):0;
342 pclCloud.at(i).normal_x = viewNormalX(0, i);
343 pclCloud.at(i).normal_y = viewNormalY(0, i);
344 pclCloud.at(i).normal_z = viewNormalZ(0, i);
349 typename PointMatcher<T>::TransformationParameters eigenMatrixToDim(
const typename PointMatcher<T>::TransformationParameters& matrix,
int dimp1)
351 typedef typename PointMatcher<T>::TransformationParameters M;
352 assert(matrix.rows() == matrix.cols());
353 assert((matrix.rows() == 3) || (matrix.rows() == 4));
354 assert((dimp1 == 3) || (dimp1 == 4));
356 if (matrix.rows() == dimp1)
359 M out(M::Identity(dimp1,dimp1));
360 out.topLeftCorner(2,2) = matrix.topLeftCorner(2,2);
361 out.topRightCorner(2,1) = matrix.topRightCorner(2,1);
366 struct KDTreeMatcherIntensity :
public PointMatcher<T>::Matcher
368 typedef PointMatcherSupport::Parametrizable Parametrizable;
369 typedef PointMatcherSupport::Parametrizable P;
370 typedef Parametrizable::Parameters Parameters;
371 typedef Parametrizable::ParameterDoc ParameterDoc;
372 typedef Parametrizable::ParametersDoc ParametersDoc;
374 typedef typename Nabo::NearestNeighbourSearch<T> NNS;
375 typedef typename NNS::SearchType NNSearchType;
377 typedef typename PointMatcher<T>::DataPoints DataPoints;
378 typedef typename PointMatcher<T>::Matcher Matcher;
379 typedef typename PointMatcher<T>::Matches Matches;
380 typedef typename PointMatcher<T>::Matrix Matrix;
384 return "This matcher matches a point from the reading to its closest neighbors in the reference.";
386 inline static const ParametersDoc availableParameters()
389 {
"knn",
"number of nearest neighbors to consider it the reference",
"1",
"1",
"2147483647", &P::Comp<unsigned>},
390 {
"epsilon",
"approximation to use for the nearest-neighbor search",
"0",
"0",
"inf", &P::Comp<T>},
391 {
"searchType",
"Nabo search type. 0: brute force, check distance to every point in the data (very slow), 1: kd-tree with linear heap, good for small knn (~up to 30) and 2: kd-tree with tree heap, good for large knn (~from 30)",
"1",
"0",
"2", &P::Comp<unsigned>},
392 {
"maxDist",
"maximum distance to consider for neighbors",
"inf",
"0",
"inf", &P::Comp<T>}
398 const NNSearchType searchType;
402 std::shared_ptr<NNS> featureNNS;
403 Matrix filteredReferenceIntensity;
406 KDTreeMatcherIntensity(
const Parameters& params = Parameters()) :
407 PointMatcher<T>::Matcher(
"KDTreeMatcherIntensity", KDTreeMatcherIntensity::availableParameters(), params),
408 knn(Parametrizable::get<int>(
"knn")),
409 epsilon(Parametrizable::get<T>(
"epsilon")),
410 searchType(NNSearchType(Parametrizable::get<int>(
"searchType"))),
411 maxDist(Parametrizable::get<T>(
"maxDist"))
413 UINFO(
"* KDTreeMatcherIntensity: initialized with knn=%d, epsilon=%f, searchType=%d and maxDist=%f", knn, epsilon, searchType, maxDist);
415 virtual ~KDTreeMatcherIntensity() {}
416 virtual void init(
const DataPoints& filteredReference)
421 filteredReferenceIntensity = filteredReference.getDescriptorCopyByName(
"intensity");
425 UWARN(
"KDTreeMatcherIntensity: knn is not over 1 (%d), intensity re-ordering will be ignored.", knn);
427 featureNNS.reset( NNS::create(filteredReference.features, filteredReference.features.rows() - 1, searchType, NNS::TOUCH_STATISTICS));
429 virtual PM::Matches findClosests(
const DP& filteredReading)
431 const int pointsCount(filteredReading.features.cols());
433 typename Matches::Dists(knn, pointsCount),
434 typename Matches::Ids(knn, pointsCount)
437 const BOOST_AUTO(filteredReadingIntensity, filteredReading.getDescriptorViewByName(
"intensity"));
439 static_assert(NNS::InvalidIndex == PM::Matches::InvalidId,
"");
440 static_assert(NNS::InvalidValue == PM::Matches::InvalidDist,
"");
441 this->visitCounter += featureNNS->knn(filteredReading.features, matches.ids, matches.dists, knn, epsilon, NNS::ALLOW_SELF_MATCH, maxDist);
445 Matches matchesOrderedByIntensity(
446 typename Matches::Dists(1, pointsCount),
447 typename Matches::Ids(1, pointsCount)
449 #pragma omp parallel for 450 for (
int i = 0; i < pointsCount; ++i)
453 for(
int k=0; k<knn && k<filteredReferenceIntensity.rows(); ++k)
455 float distIntensity = fabs(filteredReadingIntensity(0,i) - filteredReferenceIntensity(0, matches.ids.coeff(k, i)));
456 if(distIntensity < minDistance)
458 matchesOrderedByIntensity.ids.coeffRef(0, i) = matches.ids.coeff(k, i);
459 matchesOrderedByIntensity.dists.coeffRef(0, i) = matches.dists.coeff(k, i);
460 minDistance = distIntensity;
464 matches = matchesOrderedByIntensity;
475 _maxTranslation(Parameters::defaultIcpMaxTranslation()),
476 _maxRotation(Parameters::defaultIcpMaxRotation()),
477 _voxelSize(Parameters::defaultIcpVoxelSize()),
478 _downsamplingStep(Parameters::defaultIcpDownsamplingStep()),
479 _rangeMin(Parameters::defaultIcpRangeMin()),
480 _rangeMax(Parameters::defaultIcpRangeMax()),
481 _maxCorrespondenceDistance(Parameters::defaultIcpMaxCorrespondenceDistance()),
482 _maxIterations(Parameters::defaultIcpIterations()),
483 _epsilon(Parameters::defaultIcpEpsilon()),
484 _correspondenceRatio(Parameters::defaultIcpCorrespondenceRatio()),
485 _pointToPlane(Parameters::defaultIcpPointToPlane()),
486 _pointToPlaneK(Parameters::defaultIcpPointToPlaneK()),
487 _pointToPlaneRadius(Parameters::defaultIcpPointToPlaneRadius()),
488 _pointToPlaneGroundNormalsUp(Parameters::defaultIcpPointToPlaneGroundNormalsUp()),
489 _pointToPlaneMinComplexity(Parameters::defaultIcpPointToPlaneMinComplexity()),
490 _pointToPlaneLowComplexityStrategy(Parameters::defaultIcpPointToPlaneLowComplexityStrategy()),
491 _libpointmatcher(Parameters::defaultIcpPM()),
492 _libpointmatcherConfig(Parameters::defaultIcpPMConfig()),
493 _libpointmatcherKnn(Parameters::defaultIcpPMMatcherKnn()),
494 _libpointmatcherEpsilon(Parameters::defaultIcpPMMatcherEpsilon()),
495 _libpointmatcherIntensity(Parameters::defaultIcpPMMatcherIntensity()),
496 _libpointmatcherOutlierRatio(Parameters::defaultIcpPMOutlierRatio()),
497 _libpointmatcherICP(0)
504 #ifdef RTABMAP_POINTMATCHER 539 #ifndef RTABMAP_POINTMATCHER 542 UWARN(
"Parameter %s is set to true but RTAB-Map has not been built with libpointmatcher support. Setting to false.", Parameters::kIcpPM().c_str());
559 bool useDefaults =
true;
566 icp->loadFromYaml(ifs);
580 icp->readingDataPointsFilters.clear();
581 icp->readingDataPointsFilters.push_back(PM::get().DataPointsFilterRegistrar.create(
"IdentityDataPointsFilter"));
583 icp->referenceDataPointsFilters.clear();
584 icp->referenceDataPointsFilters.push_back(PM::get().DataPointsFilterRegistrar.create(
"IdentityDataPointsFilter"));
586 PM::Parameters params;
593 icp->matcher.reset(
new KDTreeMatcherIntensity<float>(params));
597 #if POINTMATCHER_VERSION_INT >= 10300 598 icp->matcher = PM::get().MatcherRegistrar.create(
"KDTreeMatcher", params);
600 icp->matcher.reset(PM::get().MatcherRegistrar.create(
"KDTreeMatcher", params));
606 icp->outlierFilters.clear();
607 icp->outlierFilters.push_back(PM::get().OutlierFilterRegistrar.create(
"TrimmedDistOutlierFilter", params));
612 icp->outlierFilters.push_back(PM::get().OutlierFilterRegistrar.create(
"SurfaceNormalOutlierFilter", params));
616 #if POINTMATCHER_VERSION_INT >= 10300 617 icp->errorMinimizer = PM::get().ErrorMinimizerRegistrar.create(
"PointToPlaneErrorMinimizer", params);
619 icp->errorMinimizer.reset(PM::get().ErrorMinimizerRegistrar.create(
"PointToPlaneErrorMinimizer", params));
625 #if POINTMATCHER_VERSION_INT >= 10300 626 icp->errorMinimizer = PM::get().ErrorMinimizerRegistrar.create(
"PointToPointErrorMinimizer");
628 icp->errorMinimizer.reset(PM::get().ErrorMinimizerRegistrar.create(
"PointToPointErrorMinimizer"));
632 icp->transformationCheckers.clear();
634 icp->transformationCheckers.push_back(PM::get().TransformationCheckerRegistrar.create(
"CounterTransformationChecker", params));
640 icp->transformationCheckers.push_back(PM::get().TransformationCheckerRegistrar.create(
"DifferentialTransformationChecker", params));
645 icp->transformationCheckers.push_back(PM::get().TransformationCheckerRegistrar.create(
"BoundTransformationChecker", params));
686 UDEBUG(
"size from=%d (format=%d, max pts=%d) to=%d (format=%d, max pts=%d)",
706 if(fromScan.
size() && toScan.
size())
709 bool hasConverged =
false;
710 float correspondencesRatio = 0.0f;
711 int correspondences = 0;
712 double variance = 1.0;
713 bool transformComputed =
false;
714 bool tooLowComplexityForPlaneToPlane =
false;
715 float secondEigenValue = 1.0f;
716 cv::Mat complexityVectors;
726 cv::Mat complexityVectorsFrom, complexityVectorsTo;
727 cv::Mat complexityValuesFrom, complexityValuesTo;
730 float complexity = fromComplexity<toComplexity?fromComplexity:toComplexity;
734 tooLowComplexityForPlaneToPlane =
true;
735 if(complexity > 0.0
f)
737 complexityVectors = fromComplexity<toComplexity?complexityVectorsFrom:complexityVectorsTo;
739 UASSERT((complexityVectors.rows == 2 && complexityVectors.cols == 2)||
740 (complexityVectors.rows == 3 && complexityVectors.cols == 3));
741 secondEigenValue = complexityValuesFrom.at<
float>(1,0)<complexityValuesTo.at<
float>(1,0)?complexityValuesFrom.at<
float>(1,0):complexityValuesTo.at<
float>(1,0);
742 UWARN(
"ICP PointToPlane ignored as structural complexity is too low (corridor-like environment): (from=%f || to=%f) < %f (%s). Second eigen value=%f. " 743 "PointToPoint is done instead, orientation is still optimized but translation will be limited to " 744 "direction of normals (%s: %s).",
747 fromComplexity<toComplexity?
"From":
"To",
748 complexityVectors.rows==2?
749 uFormat(
"n=%f,%f", complexityVectors.at<
float>(0,0), complexityVectors.at<
float>(0,1)).c_str():
751 uFormat(
"n=%f,%f,%f", complexityVectors.at<
float>(0,0), complexityVectors.at<
float>(0,1), complexityVectors.at<
float>(0,2)).c_str():
752 uFormat(
"n1=%f,%f,%f n2=%f,%f,%f", complexityVectors.at<
float>(0,0), complexityVectors.at<
float>(0,1), complexityVectors.at<
float>(0,2), complexityVectors.at<
float>(1,0), complexityVectors.at<
float>(1,1), complexityVectors.at<
float>(1,2)).c_str());
756 UWARN(
"ICP PointToPlane ignored as structural complexity cannot be computed (from=%f to=%f)!? PointToPoint is done instead.", fromComplexity, toComplexity);
760 std::cout <<
"complexityVectorsFrom = " << std::endl << complexityVectorsFrom << std::endl;
761 std::cout <<
"complexityValuesFrom = " << std::endl << complexityValuesFrom << std::endl;
762 std::cout <<
"complexityVectorsTo = " << std::endl << complexityVectorsTo << std::endl;
763 std::cout <<
"complexityValuesTo = " << std::endl << complexityValuesTo << std::endl;
774 if(fromCloudNormals->size() > 2 && toCloudNormals->size() > 2)
776 pcl::PCA<pcl::PointXYZINormal> pca;
777 pca.setInputCloud(fromCloudNormals);
778 Eigen::Vector3f valuesFrom = pca.getEigenValues();
779 pca.setInputCloud(toCloudNormals);
780 Eigen::Vector3f valuesTo = pca.getEigenValues();
781 if(valuesFrom[0]/fromCloudNormals->size() < valuesTo[0]/toCloudNormals->size())
792 pcl::PointCloud<pcl::PointXYZINormal>::Ptr fromCloudNormalsRegistered(
new pcl::PointCloud<pcl::PointXYZINormal>());
793 #ifdef RTABMAP_POINTMATCHER 797 DP data = laserScanToDP(fromScan);
801 PM::TransformationParameters T;
806 UDEBUG(
"libpointmatcher icp... (if there is a seg fault here, make sure all third party libraries are built with same Eigen version.)");
808 icpT =
Transform::fromEigen3d(Eigen::Affine3d(Eigen::Matrix4d(eigenMatrixToDim<double>(T.template cast<double>(), 4))));
811 float matchRatio = icp.errorMinimizer->getWeightedPointUsedRatio();
812 UDEBUG(
"match ratio: %f", matchRatio);
820 catch(
const std::exception &
e)
822 msg =
uFormat(
"libpointmatcher has failed: %s", e.what());
834 *fromCloudNormalsRegistered,
839 if(!icpT.
isNull() && hasConverged)
842 fromCloudNormalsRegistered,
849 transformComputed =
true;
853 int maxLaserScansFrom = fromScan.
maxPoints();
854 int maxLaserScansTo = toScan.
maxPoints();
855 if(!transformComputed)
861 if(fromCloud->size() > 2 && toCloud->size() > 2)
863 pcl::PCA<pcl::PointXYZI> pca;
864 pca.setInputCloud(fromCloud);
865 Eigen::Vector3f valuesFrom = pca.getEigenValues();
866 pca.setInputCloud(toCloud);
867 Eigen::Vector3f valuesTo = pca.getEigenValues();
868 if(valuesFrom[0]/fromCloud->size() < valuesTo[0]/toCloud->size())
876 UDEBUG(
"Computed icpStructuralDistribution %f s",timer.
ticks());
879 pcl::PointCloud<pcl::PointXYZI>::Ptr fromCloudFiltered = fromCloud;
880 pcl::PointCloud<pcl::PointXYZI>::Ptr toCloudFiltered = toCloud;
883 float pointsBeforeFiltering = (float)fromCloudFiltered->size();
885 float ratioFrom = float(fromCloudFiltered->size()) / pointsBeforeFiltering;
886 maxLaserScansFrom = int(
float(maxLaserScansFrom) * ratioFrom);
888 pointsBeforeFiltering = (float)toCloudFiltered->size();
890 float ratioTo = float(toCloudFiltered->size()) / pointsBeforeFiltering;
891 maxLaserScansTo = int(
float(maxLaserScansTo) * ratioTo);
893 UDEBUG(
"Voxel filtering time (voxel=%f m, ratioFrom=%f->%d/%d ratioTo=%f->%d/%d) = %f s",
896 (
int)fromCloudFiltered->size(),
899 (int)toCloudFiltered->size(),
904 pcl::PointCloud<pcl::PointXYZI>::Ptr fromCloudRegistered(
new pcl::PointCloud<pcl::PointXYZI>());
906 !tooLowComplexityForPlaneToPlane &&
910 pcl::PointCloud<pcl::Normal>::Ptr normalsFrom;
936 Eigen::Vector3f viewpointTo(toT.
x(), toT.
y(), toT.
z());
937 pcl::PointCloud<pcl::Normal>::Ptr normalsTo;
962 cv::Mat complexityVectorsFrom, complexityVectorsTo;
963 cv::Mat complexityValuesFrom, complexityValuesTo;
966 float complexity = fromComplexity<toComplexity?fromComplexity:toComplexity;
970 tooLowComplexityForPlaneToPlane =
true;
971 if(complexity > 0.0
f)
973 complexityVectors = fromComplexity<toComplexity?complexityVectorsFrom:complexityVectorsTo;
974 UASSERT((complexityVectors.rows == 2 && complexityVectors.cols == 2)||
975 (complexityVectors.rows == 3 && complexityVectors.cols == 3));
976 UWARN(
"ICP PointToPlane ignored as structural complexity is too low (corridor-like environment): (from=%f || to=%f) < %f (%s). " 977 "PointToPoint is done instead, orientation is still optimized but translation will be limited to " 978 "direction of normals (%s: %s).",
980 fromComplexity<toComplexity?
"From":
"To",
981 complexityVectors.rows==2?
982 uFormat(
"n=%f,%f", complexityVectors.at<
float>(0,0), complexityVectors.at<
float>(0,1)).c_str():
983 uFormat(
"n1=%f,%f,%f n2=%f,%f,%f", complexityVectors.at<
float>(0,0), complexityVectors.at<
float>(0,1), complexityVectors.at<
float>(0,2), complexityVectors.at<
float>(1,0), complexityVectors.at<
float>(1,1), complexityVectors.at<
float>(1,2)).c_str());
987 UWARN(
"ICP PointToPlane ignored as structural complexity cannot be computed (from=%f to=%f)!? PointToPoint is done instead.", fromComplexity, toComplexity);
991 std::cout <<
"complexityVectorsFrom = " << std::endl << complexityVectorsFrom << std::endl;
992 std::cout <<
"complexityValuesFrom = " << std::endl << complexityValuesFrom << std::endl;
993 std::cout <<
"complexityVectorsTo = " << std::endl << complexityVectorsTo << std::endl;
994 std::cout <<
"complexityValuesTo = " << std::endl << complexityValuesTo << std::endl;
999 pcl::PointCloud<pcl::PointXYZINormal>::Ptr fromCloudNormals(
new pcl::PointCloud<pcl::PointXYZINormal>);
1000 pcl::concatenateFields(*fromCloudFiltered, *normalsFrom, *fromCloudNormals);
1002 pcl::PointCloud<pcl::PointXYZINormal>::Ptr toCloudNormals(
new pcl::PointCloud<pcl::PointXYZINormal>);
1003 pcl::concatenateFields(*toCloudFiltered, *normalsTo, *toCloudNormals);
1005 std::vector<int> indices;
1018 Eigen::Vector3f viewpointTo(toT.
x(), toT.
y(), toT.
z()+10);
1065 UDEBUG(
"Compute normals (%d,%d) time = %f s", (
int)fromCloudNormals->size(), (int)toCloudNormals->size(), timer.
ticks());
1069 if(toCloudNormals->size() && fromCloudNormals->size())
1071 pcl::PointCloud<pcl::PointXYZINormal>::Ptr fromCloudNormalsRegistered(
new pcl::PointCloud<pcl::PointXYZINormal>());
1073 #ifdef RTABMAP_POINTMATCHER 1077 DP data = laserScanToDP(fromScan);
1081 PM::TransformationParameters T;
1086 UDEBUG(
"libpointmatcher icp... (if there is a seg fault here, make sure all third party libraries are built with same Eigen version.)");
1088 UDEBUG(
"libpointmatcher icp...done!");
1089 icpT =
Transform::fromEigen3d(Eigen::Affine3d(Eigen::Matrix4d(eigenMatrixToDim<double>(T.template cast<double>(), 4))));
1091 float matchRatio = icp.errorMinimizer->getWeightedPointUsedRatio();
1092 UDEBUG(
"match ratio: %f", matchRatio);
1097 hasConverged =
true;
1100 catch(
const std::exception &
e)
1102 msg =
uFormat(
"libpointmatcher has failed: %s", e.what());
1114 *fromCloudNormalsRegistered,
1119 if(!icpT.
isNull() && hasConverged)
1122 fromCloudNormalsRegistered,
1130 transformComputed =
true;
1134 if(!transformComputed)
1138 UWARN(
"ICP PointToPlane ignored for 2d scans with PCL registration (some crash issues). Use libpointmatcher (%s) or disable %s to avoid this warning.", Parameters::kIcpPM().c_str(), Parameters::kIcpPointToPlane().c_str());
1141 if(
_voxelSize > 0.0
f || !tooLowComplexityForPlaneToPlane)
1188 #ifdef RTABMAP_POINTMATCHER 1192 DP data = laserScanToDP(fromScan);
1196 PM::TransformationParameters T;
1201 UDEBUG(
"libpointmatcher icp... (if there is a seg fault here, make sure all third party libraries are built with same Eigen version.)");
1205 PM::ICP icpTmp =
icp;
1207 PM::Parameters params;
1209 UWARN(
"libpointmatcher icp...temporary maxDist=%s (%s=%f, %s=%f)", params[
"maxDist"].c_str(), Parameters::kIcpMaxCorrespondenceDistance().c_str(),
_maxCorrespondenceDistance, Parameters::kIcpVoxelSize().c_str(),
_voxelSize);
1212 #if POINTMATCHER_VERSION_INT >= 10300 1213 icpTmp.matcher = PM::get().MatcherRegistrar.create(
"KDTreeMatcher", params);
1215 icpTmp.matcher.reset(PM::get().MatcherRegistrar.create(
"KDTreeMatcher", params));
1218 #if POINTMATCHER_VERSION_INT >= 10300 1219 icpTmp.errorMinimizer = PM::get().ErrorMinimizerRegistrar.create(
"PointToPointErrorMinimizer");
1221 icpTmp.errorMinimizer.reset(PM::get().ErrorMinimizerRegistrar.create(
"PointToPointErrorMinimizer"));
1223 for(PM::OutlierFilters::iterator iter=icpTmp.outlierFilters.begin(); iter!=icpTmp.outlierFilters.end();)
1225 if((*iter)->className.compare(
"SurfaceNormalOutlierFilter") == 0)
1227 iter = icpTmp.outlierFilters.erase(iter);
1235 T = icpTmp(data, ref);
1241 UDEBUG(
"libpointmatcher icp...done!");
1242 icpT =
Transform::fromEigen3d(Eigen::Affine3d(Eigen::Matrix4d(eigenMatrixToDim<double>(T.template cast<double>(), 4))));
1244 float matchRatio = icp.errorMinimizer->getWeightedPointUsedRatio();
1245 UDEBUG(
"match ratio: %f", matchRatio);
1250 hasConverged =
true;
1253 catch(
const std::exception &
e)
1255 msg =
uFormat(
"libpointmatcher has failed: %s", e.what());
1267 *fromCloudRegistered,
1272 if(!icpT.
isNull() && hasConverged)
1278 msg =
uFormat(
"Rejecting transform because too low complexity (%s=0)", Parameters::kIcpPointToPlaneLowComplexityStrategy().c_str());
1286 Eigen::Vector3f v(t.
x(), t.
y(), t.
z());
1287 if(complexityVectors.cols == 2)
1290 Eigen::Vector3f n(complexityVectors.at<
float>(0,0), complexityVectors.at<
float>(0,1), 0.0f);
1292 Eigen::Vector3f vp = n*a;
1293 UWARN(
"Normals low complexity: Limiting translation from (%f,%f) to (%f,%f)",
1294 v[0], v[1], vp[0], vp[1]);
1297 else if(complexityVectors.rows == 3)
1300 Eigen::Vector3f n1(complexityVectors.at<
float>(0,0), complexityVectors.at<
float>(0,1), complexityVectors.at<
float>(0,2));
1301 Eigen::Vector3f n2(complexityVectors.at<
float>(1,0), complexityVectors.at<
float>(1,1), complexityVectors.at<
float>(1,2));
1302 float a = v.dot(n1);
1303 float b = v.dot(n2);
1304 Eigen::Vector3f vp = n1*a;
1309 UWARN(
"Normals low complexity: Limiting translation from (%f,%f,%f) to (%f,%f,%f)",
1310 v[0], v[1], v[2], vp[0], vp[1], vp[2]);
1315 UWARN(
"not supposed to be here!");
1316 v = Eigen::Vector3f(0,0,0);
1320 t =
Transform(v[0], v[1], v[2], roll, pitch, yaw);
1321 icpT = guess * t.
inverse() * guessInv;
1330 fromCloudNormalsRegistered,
1340 fromCloudRegistered,
1350 if(tooLowComplexityForPlaneToPlane)
1352 UWARN(
"Even if complexity is low , PointToPoint transformation is accepted \"as is\" (%s=2)", Parameters::kIcpPointToPlaneLowComplexityStrategy().c_str());
1355 fromCloudRegistered,
1366 if(!icpT.
isNull() && hasConverged)
1368 float ix,iy,iz, iroll,ipitch,iyaw;
1379 msg =
uFormat(
"Cannot compute transform (ICP correction too large -> %f m %f rad, limits=%f m, %f rad)",
1389 int maxLaserScans = maxLaserScansTo?maxLaserScansTo:maxLaserScansFrom;
1390 UDEBUG(
"Max scans=%d (from=%d, to=%d)", maxLaserScans, maxLaserScansFrom, maxLaserScansTo);
1394 correspondencesRatio = float(correspondences)/float(maxLaserScans);
1398 static bool warningShown =
false;
1401 UWARN(
"Maximum laser scans points not set for signature %d, correspondences ratio set relative instead of absolute! This message will only appear once.",
1403 warningShown =
true;
1405 correspondencesRatio = float(correspondences)/float(toScan.
size()>fromScan.
size()?toScan.
size():fromScan.
size());
1410 UDEBUG(
"%d->%d hasConverged=%s, variance=%f, correspondences=%d/%d (%f%%), from guess: trans=%f rot=%f",
1411 dataTo.
id(), dataFrom.
id(),
1412 hasConverged?
"true":
"false",
1415 maxLaserScans>0?maxLaserScans:(int)(toScan.
size()>fromScan.
size()?toScan.
size():fromScan.
size()),
1416 correspondencesRatio*100.0
f,
1420 if(correspondences == 0)
1422 UWARN(
"Transform is found (%s) but no correspondences has been found!? Variance is unknown!",
1427 info.
covariance = cv::Mat::eye(6,6,CV_64FC1)*variance;
1428 info.
covariance(cv::Range(3,6),cv::Range(3,6))/=10.0;
1435 msg =
uFormat(
"Cannot compute transform (cor=%d corrRatio=%f/%f maxLaserScans=%d)",
1441 transform = icpT.
inverse()*guess;
1449 msg =
uFormat(
"Cannot compute transform (converged=%s var=%f)",
1450 hasConverged?
"true":
"false", variance);
1457 msg =
"Laser scans empty ?!?";
1465 msg =
"RegistrationIcp cannot do registration with a null guess.";
1469 msg =
uFormat(
"Laser scans empty?!? (new[%d]=%d old[%d]=%d)",
GLM_FUNC_DECL T roll(detail::tquat< T, P > const &x)
static bool parse(const ParametersMap ¶meters, const std::string &key, bool &value)
int _pointToPlaneLowComplexityStrategy
cv::Point3f RTABMAP_EXP transformPoint(const cv::Point3f &pt, const Transform &transform)
const cv::Mat & data() const
pcl::PointCloud< pcl::Normal >::Ptr RTABMAP_EXP computeNormals2D(const pcl::PointCloud< pcl::PointXYZ >::Ptr &cloud, int searchK=5, float searchRadius=0.0f, const Eigen::Vector3f &viewPoint=Eigen::Vector3f(0, 0, 0))
void * _libpointmatcherICP
pcl::PointCloud< pcl::PointXYZ >::Ptr RTABMAP_EXP transformPointCloud(const pcl::PointCloud< pcl::PointXYZ >::Ptr &cloud, const Transform &transform)
void setLaserScan(const LaserScan &laserScan, bool clearPreviousData=true)
GLM_FUNC_DECL vecType< T, P > sqrt(vecType< T, P > const &x)
def init(descriptorDim, matchThreshold, iterations, cuda, model_path)
Transform RTABMAP_EXP icpPointToPlane(const pcl::PointCloud< pcl::PointNormal >::ConstPtr &cloud_source, const pcl::PointCloud< pcl::PointNormal >::ConstPtr &cloud_target, double maxCorrespondenceDistance, int maximumIterations, bool &hasConverged, pcl::PointCloud< pcl::PointNormal > &cloud_source_registered, float epsilon=0.0f, bool icp2D=false)
GLM_FUNC_DECL genType e()
float _pointToPlaneGroundNormalsUp
std::map< std::string, std::string > ParametersMap
void RTABMAP_EXP computeVarianceAndCorrespondences(const pcl::PointCloud< pcl::PointNormal >::ConstPtr &cloudA, const pcl::PointCloud< pcl::PointNormal >::ConstPtr &cloudB, double maxCorrespondenceDistance, double maxCorrespondenceAngle, double &variance, int &correspondencesOut)
int getNormalsOffset() const
Basic mathematics functions.
bool _libpointmatcherIntensity
Some conversion functions.
const LaserScan & laserScanRaw() const
GLM_FUNC_DECL T pitch(detail::tquat< T, P > const &x)
pcl::PointCloud< pcl::Normal >::Ptr RTABMAP_EXP computeFastOrganizedNormals2D(const pcl::PointCloud< pcl::PointXYZ >::Ptr &cloud, int searchK=5, float searchRadius=0.0f, const Eigen::Vector3f &viewPoint=Eigen::Vector3f(0, 0, 0))
cv::Mat RTABMAP_EXP laserScan2dFromPointCloud(const pcl::PointCloud< pcl::PointXYZ > &cloud, const Transform &transform=Transform(), bool filterNaNs=true)
bool uIsFinite(const T &value)
#define UASSERT(condition)
float _maxCorrespondenceDistance
LaserScan RTABMAP_EXP commonFiltering(const LaserScan &scan, int downsamplingStep, float rangeMin=0.0f, float rangeMax=0.0f, float voxelSize=0.0f, int normalK=0, float normalRadius=0.0f, float groundNormalsUp=0.0f)
float _correspondenceRatio
T uMax3(const T &a, const T &b, const T &c)
virtual ~RegistrationIcp()
#define UASSERT_MSG(condition, msg_str)
pcl::PointCloud< pcl::PointXYZINormal >::Ptr RTABMAP_EXP laserScanToPointCloudINormal(const LaserScan &laserScan, const Transform &transform=Transform(), float intensity=0.0f)
Transform RTABMAP_EXP icp(const pcl::PointCloud< pcl::PointXYZ >::ConstPtr &cloud_source, const pcl::PointCloud< pcl::PointXYZ >::ConstPtr &cloud_target, double maxCorrespondenceDistance, int maximumIterations, bool &hasConverged, pcl::PointCloud< pcl::PointXYZ > &cloud_source_registered, float epsilon=0.0f, bool icp2D=false)
float _pointToPlaneRadius
virtual void parseParameters(const ParametersMap ¶meters)
cv::Mat RTABMAP_EXP computeNormals(const cv::Mat &laserScan, int searchK, float searchRadius)
Transform localTransform() const
float icpStructuralDistribution
float _libpointmatcherOutlierRatio
static ULogger::Level level()
int getIntensityOffset() const
float _pointToPlaneMinComplexity
GLM_FUNC_DECL genType max(genType const &x, genType const &y)
SensorData & sensorData()
float icpStructuralComplexity
ULogger class and convenient macros.
float _libpointmatcherEpsilon
GLM_FUNC_DECL T yaw(detail::tquat< T, P > const &x)
pcl::PointCloud< pcl::PointNormal >::Ptr RTABMAP_EXP removeNaNNormalsFromPointCloud(const pcl::PointCloud< pcl::PointNormal >::Ptr &cloud)
virtual Transform computeTransformationImpl(Signature &from, Signature &to, Transform guess, RegistrationInfo &info) const
pcl::PointCloud< pcl::PointXYZI >::Ptr RTABMAP_EXP laserScanToPointCloudI(const LaserScan &laserScan, const Transform &transform=Transform(), float intensity=0.0f)
pcl::PointCloud< pcl::PointXYZ >::Ptr RTABMAP_EXP voxelize(const pcl::PointCloud< pcl::PointXYZ >::Ptr &cloud, const pcl::IndicesPtr &indices, float voxelSize)
std::string UTILITE_EXP uFormat(const char *fmt,...)
std::string UTILITE_EXP uNumber2Str(unsigned int number)
virtual void parseParameters(const ParametersMap ¶meters)
std::string _libpointmatcherConfig
float RTABMAP_EXP computeNormalsComplexity(const LaserScan &scan, const Transform &t=Transform::getIdentity(), cv::Mat *pcaEigenVectors=0, cv::Mat *pcaEigenValues=0)
bool hasIntensity() const
GLM_FUNC_DECL matType< T, P > inverse(matType< T, P > const &m)
LaserScan RTABMAP_EXP adjustNormalsToViewPoint(const LaserScan &scan, const Eigen::Vector3f &viewpoint=Eigen::Vector3f(0, 0, 0), float groundNormalsUp=0.0f)
LaserScan RTABMAP_EXP laserScanFromPointCloud(const pcl::PCLPointCloud2 &cloud, bool filterNaNs=true, bool is2D=false, const Transform &transform=Transform())