41#ifndef PCL_SAMPLE_CONSENSUS_IMPL_SAC_MODEL_TORUS_H_
42#define PCL_SAMPLE_CONSENSUS_IMPL_SAC_MODEL_TORUS_H_
45#include <pcl/sample_consensus/sac_model_torus.h>
46#include <pcl/common/concatenate.h>
49#include <unsupported/Eigen/NonLinearOptimization>
51template <
typename Po
intT,
typename Po
intNT>
56 if (samples.size() != sample_size_) {
57 PCL_ERROR(
"[pcl::SampleConsensusTorus::isSampleGood] Wrong number of samples (is "
58 "%lu, should be %lu)!\n",
64 Eigen::Vector3f n0 = Eigen::Vector3f((*normals_)[samples[0]].getNormalVector3fMap());
65 Eigen::Vector3f n1 = Eigen::Vector3f((*normals_)[samples[1]].getNormalVector3fMap());
66 Eigen::Vector3f n2 = Eigen::Vector3f((*normals_)[samples[2]].getNormalVector3fMap());
67 Eigen::Vector3f n3 = Eigen::Vector3f((*normals_)[samples[3]].getNormalVector3fMap());
70 if (std::abs((n0).cross(n1).squaredNorm()) <
71 Eigen::NumTraits<float>::dummy_precision() ||
72 std::abs((n0).cross(n2).squaredNorm()) <
73 Eigen::NumTraits<float>::dummy_precision() ||
74 std::abs((n0).cross(n3).squaredNorm()) <
75 Eigen::NumTraits<float>::dummy_precision() ||
76 std::abs((n1).cross(n2).squaredNorm()) <
77 Eigen::NumTraits<float>::dummy_precision() ||
78 std::abs((n1).cross(n3).squaredNorm()) <
79 Eigen::NumTraits<float>::dummy_precision()) {
80 PCL_ERROR(
"[pcl::SampleConsensusModelEllipse3D::isSampleGood] Sample points "
81 "normals too similar or collinear!\n");
89crossDot(
const Eigen::Vector3f& v1,
const Eigen::Vector3f& v2,
const Eigen::Vector3f& v3)
91 return v1.cross(v2).dot(v3);
95template <
typename Po
intT,
typename Po
intNT>
98 const Indices& samples, Eigen::VectorXf& model_coefficients)
const
102 if (!isSampleGood(samples)) {
103 PCL_ERROR(
"[pcl::SampleConsensusModelTorus::computeModelCoefficients] Invalid set "
104 "of samples given!\n");
109 PCL_ERROR(
"[pcl::SampleConsensusModelTorus::computeModelCoefficients] No input "
110 "dataset containing normals was given!\n");
123 const Eigen::Vector3f n0 = Eigen::Vector3f((*normals_)[samples[0]].getNormalVector3fMap());
124 const Eigen::Vector3f n1 = Eigen::Vector3f((*normals_)[samples[1]].getNormalVector3fMap());
125 const Eigen::Vector3f n2 = Eigen::Vector3f((*normals_)[samples[2]].getNormalVector3fMap());
126 const Eigen::Vector3f n3 = Eigen::Vector3f((*normals_)[samples[3]].getNormalVector3fMap());
128 const Eigen::Vector3f p0 = Eigen::Vector3f((*input_)[samples[0]].getVector3fMap());
129 const Eigen::Vector3f p1 = Eigen::Vector3f((*input_)[samples[1]].getVector3fMap());
130 const Eigen::Vector3f p2 = Eigen::Vector3f((*input_)[samples[2]].getVector3fMap());
131 const Eigen::Vector3f p3 = Eigen::Vector3f((*input_)[samples[3]].getVector3fMap());
133 const float a01 = crossDot(n0, n1, n2);
134 const float b01 = crossDot(n0, n1, n3);
135 const float a0 = crossDot(p2 - p1, n0, n2);
136 const float a1 = crossDot(p0 - p2, n1, n2);
137 const float b0 = crossDot(p3 - p1, n0, n3);
138 const float b1 = crossDot(p0 - p3, n1, n3);
139 const float a = crossDot(p0 - p2, p1 - p0, n2);
140 const float b = crossDot(p0 - p3, p1 - p0, n3);
148 const float k = -(a1 - b1 * a01 / b01) / (a0 - b0 * a01 / b01);
149 const float p = -(a - b * a01 / b01) / (a0 - b0 * a01 / b01);
157 const float _a = (b01 * k);
158 const float _b = (b01 * p + b0 * k + b1);
159 const float _c = (b0 * p + b);
161 const float eps = Eigen::NumTraits<float>::dummy_precision();
164 if ((_b * _b - 4 * _a * _c) < 0 || std::abs(a0 - b0 * a01) < eps ||
165 std::abs(b01) < eps || std::abs(_a) < eps) {
166 PCL_DEBUG(
"[pcl::SampleConsensusModelTorus::computeModelCoefficients] Can't "
167 "compute model coefficients with this method!\n");
171 const float s0 = (-_b + std::sqrt(_b * _b - 4 * _a * _c)) / (2 * _a);
172 const float s1 = (-_b - std::sqrt(_b * _b - 4 * _a * _c)) / (2 * _a);
174 float r_maj_stddev_cycle1 = std::numeric_limits<float>::max();
176 for (
float s : {s0, s1}) {
179 const float t0 = k * t1 + p;
182 Eigen::Vector3f d = ((p1 + n1 * t1) - (p0 + n0 * t0));
209 Eigen::MatrixXf A(4, 2);
210 A << d.dot(n0), 1, d.dot(n1), 1, d.dot(n2), 1, d.dot(n3), 1;
212 Eigen::Matrix<float, -1, -1>
B(4, 1);
213 B << -d.dot(p0), -d.dot(p1), -d.dot(p2), -d.dot(p3);
215 Eigen::Matrix<float, -1, -1> sol;
216#if EIGEN_VERSION_AT_LEAST(5, 0, 0)
217 sol = A.jacobiSvd<Eigen::ComputeThinU | Eigen::ComputeThinV>().solve(B);
219 sol = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B);
222 const float r_min = -sol(0);
223 const float D = sol(1);
229 const Eigen::Vector3f Pany = (p1 + n1 * t1);
231 const float lambda = (-d.dot(Pany) - D) / d.dot(d);
233 const Eigen::Vector3f centroid = Pany + d * lambda;
237 const float r_maj = std::sqrt(((p0 - r_min * n0 - centroid).squaredNorm() +
238 (p1 - r_min * n1 - centroid).squaredNorm() +
239 (p2 - r_min * n2 - centroid).squaredNorm() +
240 (p3 - r_min * n3 - centroid).squaredNorm()) /
243 const float r_maj_stddev =
244 std::sqrt((std::pow(r_maj - (p0 - r_min * n0 - centroid).norm(), 2) +
245 std::pow(r_maj - (p1 - r_min * n1 - centroid).norm(), 2) +
246 std::pow(r_maj - (p2 - r_min * n2 - centroid).norm(), 2) +
247 std::pow(r_maj - (p3 - r_min * n3 - centroid).norm(), 2)) /
250 if (r_maj_stddev < r_maj_stddev_cycle1) {
251 r_maj_stddev_cycle1 = r_maj_stddev;
257 model_coefficients.resize(model_size_);
258 model_coefficients[0] = r_maj;
259 model_coefficients[1] = r_min;
261 model_coefficients[2] = centroid[0];
262 model_coefficients[3] = centroid[1];
263 model_coefficients[4] = centroid[2];
265 model_coefficients[5] = d[0];
266 model_coefficients[6] = d[1];
267 model_coefficients[7] = d[2];
273template <
typename Po
intT,
typename Po
intNT>
276 const Eigen::VectorXf& model_coefficients, std::vector<double>& distances)
const
279 if (!isModelValid(model_coefficients)) {
284 distances.resize(indices_->size());
288 for (std::size_t i = 0; i < indices_->size(); ++i) {
289 const Eigen::Vector3f pt = (*input_)[(*indices_)[i]].getVector3fMap();
290 const Eigen::Vector3f pt_n = (*normals_)[(*indices_)[i]].getNormalVector3fMap();
292 Eigen::Vector3f torus_closest;
293 projectPointToTorus(pt, pt_n, model_coefficients, torus_closest);
295 distances[i] = (torus_closest - pt).norm();
300template <
typename Po
intT,
typename Po
intNT>
303 const Eigen::VectorXf& model_coefficients,
const double threshold,
Indices& inliers)
306 if (!isModelValid(model_coefficients)) {
310 const float squared_threshold = threshold * threshold;
312 error_sqr_dists_.clear();
313 inliers.reserve(indices_->size());
314 error_sqr_dists_.reserve(indices_->size());
316 for (std::size_t i = 0; i < indices_->size(); ++i) {
317 const Eigen::Vector3f pt = (*input_)[(*indices_)[i]].getVector3fMap();
318 const Eigen::Vector3f pt_n = (*normals_)[(*indices_)[i]].getNormalVector3fMap();
320 Eigen::Vector3f torus_closest;
321 projectPointToTorus(pt, pt_n, model_coefficients, torus_closest);
323 const float distance = (torus_closest - pt).squaredNorm();
325 if (distance < squared_threshold) {
328 inliers.push_back((*indices_)[i]);
329 error_sqr_dists_.push_back(std::sqrt(distance));
335template <
typename Po
intT,
typename Po
intNT>
338 const Eigen::VectorXf& model_coefficients,
const double threshold)
const
340 if (!isModelValid(model_coefficients))
343 const float squared_threshold = threshold * threshold;
344 std::size_t nr_p = 0;
346 for (std::size_t i = 0; i < indices_->size(); ++i) {
347 const Eigen::Vector3f pt = (*input_)[(*indices_)[i]].getVector3fMap();
348 const Eigen::Vector3f pt_n = (*normals_)[(*indices_)[i]].getNormalVector3fMap();
350 Eigen::Vector3f torus_closest;
351 projectPointToTorus(pt, pt_n, model_coefficients, torus_closest);
353 const float distance = (torus_closest - pt).squaredNorm();
355 if (distance < squared_threshold) {
363template <
typename Po
intT,
typename Po
intNT>
367 const Eigen::VectorXf& model_coefficients,
368 Eigen::VectorXf& optimized_coefficients)
const
371 optimized_coefficients = model_coefficients;
374 if (!isModelValid(model_coefficients)) {
375 PCL_ERROR(
"[pcl::SampleConsensusModelTorus::optimizeModelCoefficients] Given model "
381 if (inliers.size() <= sample_size_) {
382 PCL_ERROR(
"[pcl::SampleConsensusModelTorus::optimizeModelCoefficients] Not enough "
383 "inliers to refine/optimize the model's coefficients (%lu)! Returning "
384 "the same coefficients.\n",
389 OptimizationFunctor functor(
this, inliers);
390 Eigen::NumericalDiff<OptimizationFunctor> num_diff(functor);
391 Eigen::LevenbergMarquardt<Eigen::NumericalDiff<OptimizationFunctor>,
double> lm(
394 Eigen::VectorXd coeff = model_coefficients.cast<
double>();
395 int info = lm.minimize(coeff);
399 "[pcl::SampleConsensusModelTorus::optimizeModelCoefficients] Optimizer returned"
400 "with error (%i)! Returning ",
406 coeff.tail<3>().normalize();
407 optimized_coefficients = coeff.cast<
float>();
411template <
typename Po
intT,
typename Po
intNT>
414 const Eigen::Vector3f& p_in,
415 const Eigen::Vector3f& p_n,
416 const Eigen::VectorXf& model_coefficients,
417 Eigen::Vector3f& pt_out)
const
421 const float& R = model_coefficients[0];
422 const float& r = model_coefficients[1];
424 const float& x0 = model_coefficients[2];
425 const float& y0 = model_coefficients[3];
426 const float& z0 = model_coefficients[4];
428 const float& nx = model_coefficients[5];
429 const float& ny = model_coefficients[6];
430 const float& nz = model_coefficients[7];
433 Eigen::Vector3f n{nx, ny, nz};
436 Eigen::Vector3f pt0{x0, y0, z0};
439 const float D = -n.dot(pt0);
449 Eigen::Vector3f pt_proj;
453 if (std::abs(n.dot(p_n)) > Eigen::NumTraits<float>::dummy_precision()) {
454 float lambda = (-D - n.dot(p_in)) / n.dot(p_n);
455 pt_proj = p_in + lambda * p_n;
462 const Eigen::Vector3f circle_closest = (pt_proj - pt0).normalized() * R + pt0;
466 pt_out = (p_in - circle_closest).normalized() * r + circle_closest;
470template <
typename Po
intT,
typename Po
intNT>
474 const Eigen::VectorXf& model_coefficients,
475 PointCloud& projected_points,
476 bool copy_data_fields)
const
479 if (!isModelValid(model_coefficients)) {
481 "[pcl::SampleConsensusModelCylinder::projectPoints] Given model is invalid!\n");
486 if (copy_data_fields) {
488 projected_points.resize(input_->size());
489 projected_points.width = input_->width;
490 projected_points.height = input_->height;
492 using FieldList =
typename pcl::traits::fieldList<PointT>::type;
494 for (std::size_t i = 0; i < input_->size(); ++i)
496 pcl::for_each_type<FieldList>(
500 for (
const auto& inlier : inliers) {
502 const Eigen::Vector3f pt_n = (*normals_)[inlier].getNormalVector3fMap();
504 (*input_)[inlier].getVector3fMap(), pt_n, model_coefficients, q);
505 projected_points[inlier].getVector3fMap() = q;
510 projected_points.resize(inliers.size());
511 projected_points.width = inliers.size();
512 projected_points.height = 1;
514 using FieldList =
typename pcl::traits::fieldList<PointT>::type;
516 for (std::size_t i = 0; i < inliers.size(); ++i) {
519 (*input_)[inliers[i]], projected_points[i]));
522 for (
const auto& inlier : inliers) {
524 const Eigen::Vector3f pt_n = (*normals_)[inlier].getNormalVector3fMap();
526 (*input_)[inlier].getVector3fMap(), pt_n, model_coefficients, q);
527 projected_points[inlier].getVector3fMap() = q;
533template <
typename Po
intT,
typename Po
intNT>
536 const std::set<index_t>& indices,
537 const Eigen::VectorXf& model_coefficients,
538 const double threshold)
const
541 for (
const auto& index : indices) {
542 const Eigen::Vector3f pt_n = (*normals_)[index].getNormalVector3fMap();
543 Eigen::Vector3f torus_closest;
544 projectPointToTorus((*input_)[index].getVector3fMap(), pt_n, model_coefficients, torus_closest);
546 if (((*input_)[index].getVector3fMap() - torus_closest).squaredNorm() > threshold * threshold)
553template <
typename Po
intT,
typename Po
intNT>
556 const Eigen::VectorXf& model_coefficients)
const
561 if (radius_min_ != std::numeric_limits<double>::lowest() &&
562 (model_coefficients[0] < radius_min_ || model_coefficients[1] < radius_min_)) {
564 "[pcl::SampleConsensusModelTorus::isModelValid] Major radius OR minor radius "
565 "of torus is/are too small: should be larger than %g, but are {%g, %g}.\n",
567 model_coefficients[0],
568 model_coefficients[1]);
571 if (radius_max_ != std::numeric_limits<double>::max() &&
572 (model_coefficients[0] > radius_max_ || model_coefficients[1] > radius_max_)) {
574 "[pcl::SampleConsensusModelTorus::isModelValid] Major radius OR minor radius "
575 "of torus is/are too big: should be smaller than %g, but are {%g, %g}.\n",
577 model_coefficients[0],
578 model_coefficients[1]);
584#define PCL_INSTANTIATE_SampleConsensusModelTorus(PointT, PointNT) \
585 template class PCL_EXPORTS pcl::SampleConsensusModelTorus<PointT, PointNT>;
SampleConsensusModel represents the base model class.
void projectPointToTorus(const Eigen::Vector3f &pt, const Eigen::Vector3f &pt_n, const Eigen::VectorXf &model_coefficients, Eigen::Vector3f &pt_proj) const
Project a point onto a torus given by its model coefficients (radii, torus_center_point,...
void selectWithinDistance(const Eigen::VectorXf &model_coefficients, const double threshold, Indices &inliers) override
Select all the points which respect the given model coefficients as inliers.
bool doSamplesVerifyModel(const std::set< index_t > &indices, const Eigen::VectorXf &model_coefficients, const double threshold) const override
Verify whether a subset of indices verifies the given torus model coefficients.
bool isModelValid(const Eigen::VectorXf &model_coefficients) const override
Check whether a model is valid given the user constraints.
bool computeModelCoefficients(const Indices &samples, Eigen::VectorXf &model_coefficients) const override
Check whether the given index samples can form a valid torus model, compute the model coefficients fr...
void optimizeModelCoefficients(const Indices &inliers, const Eigen::VectorXf &model_coefficients, Eigen::VectorXf &optimized_coefficients) const override
Recompute the torus coefficients using the given inlier set and return them to the user.
std::size_t countWithinDistance(const Eigen::VectorXf &model_coefficients, const double threshold) const override
Count all the points which respect the given model coefficients as inliers.
bool isSampleGood(const Indices &samples) const override
Check if a sample of indices results in a good sample of points indices.
void getDistancesToModel(const Eigen::VectorXf &model_coefficients, std::vector< double > &distances) const override
Compute all distances from the cloud data to a given torus model.
void projectPoints(const Indices &inliers, const Eigen::VectorXf &model_coefficients, PointCloud &projected_points, bool copy_data_fields=true) const override
Create a new point cloud with inliers projected onto the torus model.
IndicesAllocator<> Indices
Type used for indices in PCL.
Helper functor structure for concatenate.