41 #ifndef PCL_SAMPLE_CONSENSUS_IMPL_MLESAC_H_
42 #define PCL_SAMPLE_CONSENSUS_IMPL_MLESAC_H_
44 #include <pcl/sample_consensus/mlesac.h>
49 template <
typename Po
intT>
bool
53 if (threshold_ == std::numeric_limits<double>::max())
55 PCL_ERROR (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] No threshold set!\n");
60 double d_best_penalty = std::numeric_limits<double>::max();
64 Eigen::VectorXf model_coefficients (sac_model_->getModelSize ());
65 std::vector<double> distances;
68 sigma_ = computeMedianAbsoluteDeviation (sac_model_->getInputCloud (), sac_model_->getIndices (), threshold_);
69 const double dist_scaling_factor = -1.0 / (2.0 * sigma_ * sigma_);
70 const double normalization_factor = 1.0 / (sqrt (2 *
M_PI) * sigma_);
71 if (debug_verbosity_level > 1)
72 PCL_DEBUG (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] Estimated sigma value: %f.\n", sigma_);
75 Eigen::Vector4f min_pt, max_pt;
76 getMinMax (sac_model_->getInputCloud (), sac_model_->getIndices (), min_pt, max_pt);
78 double v = sqrt (max_pt.dot (max_pt));
80 int n_inliers_count = 0;
81 std::size_t indices_size;
82 unsigned skipped_count = 0;
84 const unsigned max_skip = max_iterations_ * 10;
87 while (iterations_ < k && skipped_count < max_skip)
90 sac_model_->getSamples (iterations_, selection);
92 if (selection.empty ())
break;
95 if (!sac_model_->computeModelCoefficients (selection, model_coefficients))
103 sac_model_->getDistancesToModel (model_coefficients, distances);
105 if (distances.empty ())
115 double p_outlier_prob = 0;
117 indices_size = sac_model_->getIndices ()->size ();
118 std::vector<double> p_inlier_prob (indices_size);
119 for (
int j = 0; j < iterations_EM_; ++j)
121 const double weighted_normalization_factor = gamma * normalization_factor;
123 for (std::size_t i = 0; i < indices_size; ++i)
124 p_inlier_prob[i] = weighted_normalization_factor * std::exp ( dist_scaling_factor * distances[i] * distances[i] );
127 p_outlier_prob = (1 - gamma) / v;
130 for (std::size_t i = 0; i < indices_size; ++i)
131 gamma += p_inlier_prob [i] / (p_inlier_prob[i] + p_outlier_prob);
132 gamma /=
static_cast<double>(sac_model_->getIndices ()->size ());
136 double d_cur_penalty = 0;
137 for (std::size_t i = 0; i < indices_size; ++i)
138 d_cur_penalty += std::log (p_inlier_prob[i] + p_outlier_prob);
139 d_cur_penalty = - d_cur_penalty;
142 if (d_cur_penalty < d_best_penalty)
144 d_best_penalty = d_cur_penalty;
148 model_coefficients_ = model_coefficients;
152 for (
const double &
distance : distances)
157 double w =
static_cast<double> (n_inliers_count) /
static_cast<double> (sac_model_->getIndices ()->size ());
158 double p_no_outliers = 1 - std::pow (w,
static_cast<double> (selection.size ()));
159 p_no_outliers = (std::max) (std::numeric_limits<double>::epsilon (), p_no_outliers);
160 p_no_outliers = (std::min) (1 - std::numeric_limits<double>::epsilon (), p_no_outliers);
161 k = std::log (1 - probability_) / std::log (p_no_outliers);
165 if (debug_verbosity_level > 1)
166 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);
167 if (iterations_ > max_iterations_)
169 if (debug_verbosity_level > 0)
170 PCL_DEBUG (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] MLESAC reached the maximum number of trials.\n");
177 if (debug_verbosity_level > 0)
178 PCL_DEBUG (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] Unable to find a solution!\n");
183 sac_model_->getDistancesToModel (model_coefficients_, distances);
184 Indices &indices = *sac_model_->getIndices ();
185 if (distances.size () != indices.size ())
187 PCL_ERROR (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] Estimated distances (%lu) differs than the normal of indices (%lu).\n", distances.size (), indices.size ());
191 inliers_.resize (distances.size ());
194 for (std::size_t i = 0; i < distances.size (); ++i)
195 if (distances[i] <= 2 * sigma_)
196 inliers_[n_inliers_count++] = indices[i];
199 inliers_.resize (n_inliers_count);
201 if (debug_verbosity_level > 0)
202 PCL_DEBUG (
"[pcl::MaximumLikelihoodSampleConsensus::computeModel] Model: %lu size, %d inliers.\n", model_.size (), n_inliers_count);
208 template <
typename Po
intT>
double
210 const PointCloudConstPtr &cloud,
214 std::vector<double> distances (indices->size ());
216 Eigen::Vector4f median;
220 for (std::size_t i = 0; i < indices->size (); ++i)
223 Eigen::Vector4f ptdiff = pt - median;
225 distances[i] = ptdiff.dot (ptdiff);
228 const double result =
pcl::computeMedian (distances.begin (), distances.end (),
static_cast<double(*)(
double)
>(std::sqrt));
229 return (sigma * result);
233 template <
typename Po
intT>
void
235 const PointCloudConstPtr &cloud,
237 Eigen::Vector4f &min_p,
238 Eigen::Vector4f &max_p)
const
240 min_p.setConstant (FLT_MAX);
241 max_p.setConstant (-FLT_MAX);
242 min_p[3] = max_p[3] = 0;
244 for (std::size_t i = 0; i < indices->size (); ++i)
246 if ((*cloud)[(*indices)[i]].x < min_p[0]) min_p[0] = (*cloud)[(*indices)[i]].x;
247 if ((*cloud)[(*indices)[i]].y < min_p[1]) min_p[1] = (*cloud)[(*indices)[i]].y;
248 if ((*cloud)[(*indices)[i]].z < min_p[2]) min_p[2] = (*cloud)[(*indices)[i]].z;
250 if ((*cloud)[(*indices)[i]].x > max_p[0]) max_p[0] = (*cloud)[(*indices)[i]].x;
251 if ((*cloud)[(*indices)[i]].y > max_p[1]) max_p[1] = (*cloud)[(*indices)[i]].y;
252 if ((*cloud)[(*indices)[i]].z > max_p[2]) max_p[2] = (*cloud)[(*indices)[i]].z;
257 template <
typename Po
intT>
void
259 const PointCloudConstPtr &cloud,
261 Eigen::Vector4f &median)
const
264 std::vector<float> x (indices->size ());
265 std::vector<float> y (indices->size ());
266 std::vector<float> z (indices->size ());
267 for (std::size_t i = 0; i < indices->size (); ++i)
269 x[i] = (*cloud)[(*indices)[i]].x;
270 y[i] = (*cloud)[(*indices)[i]].y;
271 z[i] = (*cloud)[(*indices)[i]].z;
280 #define PCL_INSTANTIATE_MaximumLikelihoodSampleConsensus(T) template class PCL_EXPORTS pcl::MaximumLikelihoodSampleConsensus<T>;
282 #endif // PCL_SAMPLE_CONSENSUS_IMPL_MLESAC_H_