41 #ifndef PCL_SAMPLE_CONSENSUS_IMPL_MLESAC_H_
42 #define PCL_SAMPLE_CONSENSUS_IMPL_MLESAC_H_
44 #include <pcl/sample_consensus/mlesac.h>
48 template <
typename Po
intT>
bool
52 if (threshold_ == std::numeric_limits<double>::max())
54 PCL_ERROR (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] No threshold set!\n");
59 double d_best_penalty = std::numeric_limits<double>::max();
63 Eigen::VectorXf model_coefficients;
64 std::vector<double> distances;
67 sigma_ = computeMedianAbsoluteDeviation (sac_model_->getInputCloud (), sac_model_->getIndices (), threshold_);
68 if (debug_verbosity_level > 1)
69 PCL_DEBUG (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] Estimated sigma value: %f.\n", sigma_);
72 Eigen::Vector4f min_pt, max_pt;
73 getMinMax (sac_model_->getInputCloud (), sac_model_->getIndices (), min_pt, max_pt);
75 double v = sqrt (max_pt.dot (max_pt));
77 int n_inliers_count = 0;
78 std::size_t indices_size;
79 unsigned skipped_count = 0;
81 const unsigned max_skip = max_iterations_ * 10;
84 while (iterations_ < k && skipped_count < max_skip)
87 sac_model_->getSamples (iterations_, selection);
89 if (selection.empty ())
break;
92 if (!sac_model_->computeModelCoefficients (selection, model_coefficients))
100 sac_model_->getDistancesToModel (model_coefficients, distances);
102 if (distances.empty ())
112 double p_outlier_prob = 0;
114 indices_size = sac_model_->getIndices ()->size ();
115 std::vector<double> p_inlier_prob (indices_size);
116 for (
int j = 0; j < iterations_EM_; ++j)
119 for (std::size_t i = 0; i < indices_size; ++i)
120 p_inlier_prob[i] = gamma * std::exp (- (distances[i] * distances[i] ) / 2 * (sigma_ * sigma_) ) /
121 (sqrt (2 *
M_PI) * sigma_);
124 p_outlier_prob = (1 - gamma) / v;
127 for (std::size_t i = 0; i < indices_size; ++i)
128 gamma += p_inlier_prob [i] / (p_inlier_prob[i] + p_outlier_prob);
129 gamma /=
static_cast<double>(sac_model_->getIndices ()->size ());
133 double d_cur_penalty = 0;
134 for (std::size_t i = 0; i < indices_size; ++i)
135 d_cur_penalty += std::log (p_inlier_prob[i] + p_outlier_prob);
136 d_cur_penalty = - d_cur_penalty;
139 if (d_cur_penalty < d_best_penalty)
141 d_best_penalty = d_cur_penalty;
145 model_coefficients_ = model_coefficients;
149 for (
const double &
distance : distances)
154 double w =
static_cast<double> (n_inliers_count) /
static_cast<double> (sac_model_->getIndices ()->size ());
155 double p_no_outliers = 1 - std::pow (w,
static_cast<double> (selection.size ()));
156 p_no_outliers = (std::max) (std::numeric_limits<double>::epsilon (), p_no_outliers);
157 p_no_outliers = (std::min) (1 - std::numeric_limits<double>::epsilon (), p_no_outliers);
158 k = std::log (1 - probability_) / std::log (p_no_outliers);
162 if (debug_verbosity_level > 1)
163 PCL_DEBUG (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] Trial %d out of %d. Best penalty is %f.\n", iterations_,
static_cast<int> (std::ceil (k)), d_best_penalty);
164 if (iterations_ > max_iterations_)
166 if (debug_verbosity_level > 0)
167 PCL_DEBUG (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] MLESAC reached the maximum number of trials.\n");
174 if (debug_verbosity_level > 0)
175 PCL_DEBUG (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] Unable to find a solution!\n");
180 sac_model_->getDistancesToModel (model_coefficients_, distances);
181 Indices &indices = *sac_model_->getIndices ();
182 if (distances.size () != indices.size ())
184 PCL_ERROR (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] Estimated distances (%lu) differs than the normal of indices (%lu).\n", distances.size (), indices.size ());
188 inliers_.resize (distances.size ());
191 for (std::size_t i = 0; i < distances.size (); ++i)
192 if (distances[i] <= 2 * sigma_)
193 inliers_[n_inliers_count++] = indices[i];
196 inliers_.resize (n_inliers_count);
198 if (debug_verbosity_level > 0)
199 PCL_DEBUG (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] Model: %lu size, %d inliers.\n", model_.size (), n_inliers_count);
205 template <
typename Po
intT>
double
207 const PointCloudConstPtr &cloud,
211 std::vector<double> distances (indices->size ());
213 Eigen::Vector4f median;
215 computeMedian (cloud, indices, median);
217 for (std::size_t i = 0; i < indices->size (); ++i)
220 Eigen::Vector4f ptdiff = pt - median;
222 distances[i] = ptdiff.dot (ptdiff);
225 std::sort (distances.begin (), distances.end ());
228 std::size_t mid = indices->size () / 2;
230 if (indices->size () % 2 == 0)
231 result = (sqrt (distances[mid-1]) + sqrt (distances[mid])) / 2;
233 result = sqrt (distances[mid]);
234 return (sigma * result);
238 template <
typename Po
intT>
void
240 const PointCloudConstPtr &cloud,
242 Eigen::Vector4f &min_p,
243 Eigen::Vector4f &max_p)
const
245 min_p.setConstant (FLT_MAX);
246 max_p.setConstant (-FLT_MAX);
247 min_p[3] = max_p[3] = 0;
249 for (std::size_t i = 0; i < indices->size (); ++i)
251 if ((*cloud)[(*indices)[i]].x < min_p[0]) min_p[0] = (*cloud)[(*indices)[i]].x;
252 if ((*cloud)[(*indices)[i]].y < min_p[1]) min_p[1] = (*cloud)[(*indices)[i]].y;
253 if ((*cloud)[(*indices)[i]].z < min_p[2]) min_p[2] = (*cloud)[(*indices)[i]].z;
255 if ((*cloud)[(*indices)[i]].x > max_p[0]) max_p[0] = (*cloud)[(*indices)[i]].x;
256 if ((*cloud)[(*indices)[i]].y > max_p[1]) max_p[1] = (*cloud)[(*indices)[i]].y;
257 if ((*cloud)[(*indices)[i]].z > max_p[2]) max_p[2] = (*cloud)[(*indices)[i]].z;
262 template <
typename Po
intT>
void
264 const PointCloudConstPtr &cloud,
266 Eigen::Vector4f &median)
const
269 std::vector<float> x (indices->size ());
270 std::vector<float> y (indices->size ());
271 std::vector<float> z (indices->size ());
272 for (std::size_t i = 0; i < indices->size (); ++i)
274 x[i] = (*cloud)[(*indices)[i]].x;
275 y[i] = (*cloud)[(*indices)[i]].y;
276 z[i] = (*cloud)[(*indices)[i]].z;
278 std::sort (x.begin (), x.end ());
279 std::sort (y.begin (), y.end ());
280 std::sort (z.begin (), z.end ());
282 std::size_t mid = indices->size () / 2;
283 if (indices->size () % 2 == 0)
285 median[0] = (x[mid-1] + x[mid]) / 2;
286 median[1] = (y[mid-1] + y[mid]) / 2;
287 median[2] = (z[mid-1] + z[mid]) / 2;
298 #define PCL_INSTANTIATE_MaximumLikelihoodSampleConsensus(T) template class PCL_EXPORTS pcl::MaximumLikelihoodSampleConsensus<T>;
void computeMedian(const PointCloudConstPtr &cloud, const IndicesPtr &indices, Eigen::Vector4f &median) const
Compute the median value of a 3D point cloud using a given set point indices and return it as a Point...
void getMinMax(const PointCloudConstPtr &cloud, const IndicesPtr &indices, Eigen::Vector4f &min_p, Eigen::Vector4f &max_p) const
Determine the minimum and maximum 3D bounding box coordinates for a given set of points.
bool computeModel(int debug_verbosity_level=0) override
Compute the actual model and find the inliers.
double computeMedianAbsoluteDeviation(const PointCloudConstPtr &cloud, const IndicesPtr &indices, double sigma) const
Compute the median absolute deviation:
Defines all the PCL implemented PointT point type structures.
void getMinMax(const PointT &histogram, int len, float &min_p, float &max_p)
Get the minimum and maximum values on a point histogram.
float distance(const PointT &p1, const PointT &p2)
const Eigen::Map< const Eigen::Vector4f, Eigen::Aligned > Vector4fMapConst
IndicesAllocator<> Indices
Type used for indices in PCL.
shared_ptr< Indices > IndicesPtr