38 #ifndef PCL_SAMPLE_CONSENSUS_IMPL_SAC_MODEL_CONE_H_
39 #define PCL_SAMPLE_CONSENSUS_IMPL_SAC_MODEL_CONE_H_
43 #include <unsupported/Eigen/NonLinearOptimization>
46 template <
typename Po
intT,
typename Po
intNT>
bool
53 template <
typename Po
intT,
typename Po
intNT>
bool
55 const std::vector<int> &samples, Eigen::VectorXf &model_coefficients)
58 if (samples.size () != 3)
60 PCL_ERROR (
"[pcl::SampleConsensusModelCone::computeModelCoefficients] Invalid set of samples given (%zu)!\n", samples.size ());
66 PCL_ERROR (
"[pcl::SampleConsensusModelCone::computeModelCoefficients] No input dataset containing normals was given!\n");
70 Eigen::Vector4f p1 (input_->points[samples[0]].x, input_->points[samples[0]].y, input_->points[samples[0]].z, 0);
71 Eigen::Vector4f p2 (input_->points[samples[1]].x, input_->points[samples[1]].y, input_->points[samples[1]].z, 0);
72 Eigen::Vector4f p3 (input_->points[samples[2]].x, input_->points[samples[2]].y, input_->points[samples[2]].z, 0);
74 Eigen::Vector4f n1 (normals_->points[samples[0]].normal[0], normals_->points[samples[0]].normal[1], normals_->points[samples[0]].normal[2], 0);
75 Eigen::Vector4f n2 (normals_->points[samples[1]].normal[0], normals_->points[samples[1]].normal[1], normals_->points[samples[1]].normal[2], 0);
76 Eigen::Vector4f n3 (normals_->points[samples[2]].normal[0], normals_->points[samples[2]].normal[1], normals_->points[samples[2]].normal[2], 0);
79 Eigen::Vector4f ortho12 = n1.cross3(n2);
80 Eigen::Vector4f ortho23 = n2.cross3(n3);
81 Eigen::Vector4f ortho31 = n3.cross3(n1);
83 float denominator = n1.dot(ortho23);
85 float d1 = p1.dot (n1);
86 float d2 = p2.dot (n2);
87 float d3 = p3.dot (n3);
89 Eigen::Vector4f apex = (d1 * ortho23 + d2 * ortho31 + d3 * ortho12) / denominator;
92 Eigen::Vector4f ap1 = p1 - apex;
93 Eigen::Vector4f ap2 = p2 - apex;
94 Eigen::Vector4f ap3 = p3 - apex;
96 Eigen::Vector4f np1 = apex + (ap1/ap1.norm ());
97 Eigen::Vector4f np2 = apex + (ap2/ap2.norm ());
98 Eigen::Vector4f np3 = apex + (ap3/ap3.norm ());
100 Eigen::Vector4f np1np2 = np2 - np1;
101 Eigen::Vector4f np1np3 = np3 - np1;
103 Eigen::Vector4f axis_dir = np1np2.cross3 (np1np3);
104 axis_dir.normalize ();
112 float opening_angle = ( acosf (ap1.dot (axis_dir)) + acosf (ap2.dot (axis_dir)) + acosf (ap3.dot (axis_dir)) ) / 3.0f;
114 model_coefficients.resize (7);
116 model_coefficients[0] = apex[0];
117 model_coefficients[1] = apex[1];
118 model_coefficients[2] = apex[2];
120 model_coefficients[3] = axis_dir[0];
121 model_coefficients[4] = axis_dir[1];
122 model_coefficients[5] = axis_dir[2];
124 model_coefficients[6] = opening_angle;
126 if (model_coefficients[6] != -std::numeric_limits<double>::max() && model_coefficients[6] < min_angle_)
128 if (model_coefficients[6] != std::numeric_limits<double>::max() && model_coefficients[6] > max_angle_)
135 template <
typename Po
intT,
typename Po
intNT>
void
137 const Eigen::VectorXf &model_coefficients, std::vector<double> &distances)
140 if (!isModelValid (model_coefficients))
146 distances.resize (indices_->size ());
148 Eigen::Vector4f apex (model_coefficients[0], model_coefficients[1], model_coefficients[2], 0);
149 Eigen::Vector4f axis_dir (model_coefficients[3], model_coefficients[4], model_coefficients[5], 0);
150 float opening_angle = model_coefficients[6];
152 float apexdotdir = apex.dot (axis_dir);
153 float dirdotdir = 1.0f / axis_dir.dot (axis_dir);
155 for (
size_t i = 0; i < indices_->size (); ++i)
157 Eigen::Vector4f pt (input_->points[(*indices_)[i]].x, input_->points[(*indices_)[i]].y, input_->points[(*indices_)[i]].z, 0);
158 Eigen::Vector4f n (normals_->points[(*indices_)[i]].normal[0], normals_->points[(*indices_)[i]].normal[1], normals_->points[(*indices_)[i]].normal[2], 0);
161 float k = (pt.dot (axis_dir) - apexdotdir) * dirdotdir;
162 Eigen::Vector4f pt_proj = apex + k * axis_dir;
163 Eigen::Vector4f dir = pt - pt_proj;
167 Eigen::Vector4f height = apex - pt_proj;
168 float actual_cone_radius = tanf (opening_angle) * height.norm ();
172 Eigen::Vector4f cone_normal = sinf (opening_angle) * height + cosf (opening_angle) * dir;
176 double d_euclid = fabs (pointToAxisDistance (pt, model_coefficients) - actual_cone_radius);
179 double d_normal = fabs (
getAngle3D (n, cone_normal));
180 d_normal = (std::min) (d_normal, M_PI - d_normal);
182 distances[i] = fabs (normal_distance_weight_ * d_normal + (1 - normal_distance_weight_) * d_euclid);
187 template <
typename Po
intT,
typename Po
intNT>
void
189 const Eigen::VectorXf &model_coefficients,
const double threshold, std::vector<int> &inliers)
192 if (!isModelValid (model_coefficients))
199 inliers.resize (indices_->size ());
201 Eigen::Vector4f apex (model_coefficients[0], model_coefficients[1], model_coefficients[2], 0);
202 Eigen::Vector4f axis_dir (model_coefficients[3], model_coefficients[4], model_coefficients[5], 0);
203 float opening_angle = model_coefficients[6];
205 float apexdotdir = apex.dot (axis_dir);
206 float dirdotdir = 1.0f / axis_dir.dot (axis_dir);
208 for (
size_t i = 0; i < indices_->size (); ++i)
210 Eigen::Vector4f pt (input_->points[(*indices_)[i]].x, input_->points[(*indices_)[i]].y, input_->points[(*indices_)[i]].z, 0);
211 Eigen::Vector4f n (normals_->points[(*indices_)[i]].normal[0], normals_->points[(*indices_)[i]].normal[1], normals_->points[(*indices_)[i]].normal[2], 0);
214 float k = (pt.dot (axis_dir) - apexdotdir) * dirdotdir;
215 Eigen::Vector4f pt_proj = apex + k * axis_dir;
218 Eigen::Vector4f pp_pt_dir = pt - pt_proj;
219 pp_pt_dir.normalize ();
222 Eigen::Vector4f height = apex - pt_proj;
223 double actual_cone_radius = tan(opening_angle) * height.norm ();
227 Eigen::Vector4f cone_normal = sinf (opening_angle) * height + cosf (opening_angle) * pp_pt_dir;
231 double d_euclid = fabs (pointToAxisDistance (pt, model_coefficients) - actual_cone_radius);
234 double d_normal = fabs (
getAngle3D (n, cone_normal));
235 d_normal = (std::min) (d_normal, M_PI - d_normal);
237 if (fabs (normal_distance_weight_ * d_normal + (1 - normal_distance_weight_) * d_euclid) < threshold)
240 inliers[nr_p] = (*indices_)[i];
244 inliers.resize (nr_p);
248 template <
typename Po
intT,
typename Po
intNT>
int
250 const Eigen::VectorXf &model_coefficients,
const double threshold)
254 if (!isModelValid (model_coefficients))
259 Eigen::Vector4f apex (model_coefficients[0], model_coefficients[1], model_coefficients[2], 0);
260 Eigen::Vector4f axis_dir (model_coefficients[3], model_coefficients[4], model_coefficients[5], 0);
261 float opening_angle = model_coefficients[6];
263 float apexdotdir = apex.dot (axis_dir);
264 float dirdotdir = 1.0f / axis_dir.dot (axis_dir);
266 for (
size_t i = 0; i < indices_->size (); ++i)
268 Eigen::Vector4f pt (input_->points[(*indices_)[i]].x, input_->points[(*indices_)[i]].y, input_->points[(*indices_)[i]].z, 0);
269 Eigen::Vector4f n (normals_->points[(*indices_)[i]].normal[0], normals_->points[(*indices_)[i]].normal[1], normals_->points[(*indices_)[i]].normal[2], 0);
272 float k = (pt.dot (axis_dir) - apexdotdir) * dirdotdir;
273 Eigen::Vector4f pt_proj = apex + k * axis_dir;
276 Eigen::Vector4f pp_pt_dir = pt - pt_proj;
277 pp_pt_dir.normalize ();
280 Eigen::Vector4f height = apex - pt_proj;
281 double actual_cone_radius = tan(opening_angle) * height.norm ();
285 Eigen::Vector4f cone_normal = sinf (opening_angle) * height + cosf (opening_angle) * pp_pt_dir;
289 double d_euclid = fabs (pointToAxisDistance (pt, model_coefficients) - actual_cone_radius);
292 double d_normal = fabs (
getAngle3D (n, cone_normal));
293 d_normal = (std::min) (d_normal, M_PI - d_normal);
295 if (fabs (normal_distance_weight_ * d_normal + (1 - normal_distance_weight_) * d_euclid) < threshold)
302 template <
typename Po
intT,
typename Po
intNT>
void
304 const std::vector<int> &inliers,
const Eigen::VectorXf &model_coefficients, Eigen::VectorXf &optimized_coefficients)
306 optimized_coefficients = model_coefficients;
309 if (model_coefficients.size () != 7)
311 PCL_ERROR (
"[pcl::SampleConsensusModelCone::optimizeModelCoefficients] Invalid number of model coefficients given (%zu)!\n", model_coefficients.size ());
315 if (inliers.empty ())
317 PCL_DEBUG (
"[pcl::SampleConsensusModelCone:optimizeModelCoefficients] Inliers vector empty! Returning the same coefficients.\n");
321 tmp_inliers_ = &inliers;
323 OptimizationFunctor functor (static_cast<int> (inliers.size ()),
this);
324 Eigen::NumericalDiff<OptimizationFunctor > num_diff (functor);
325 Eigen::LevenbergMarquardt<Eigen::NumericalDiff<OptimizationFunctor>,
float> lm (num_diff);
326 int info = lm.minimize (optimized_coefficients);
329 PCL_DEBUG (
"[pcl::SampleConsensusModelCone::optimizeModelCoefficients] LM solver finished with exit code %i, having a residual norm of %g. \nInitial solution: %g %g %g %g %g %g %g \nFinal solution: %g %g %g %g %g %g %g\n",
330 info, lm.fvec.norm (), model_coefficients[0], model_coefficients[1], model_coefficients[2], model_coefficients[3],
331 model_coefficients[4], model_coefficients[5], model_coefficients[6], optimized_coefficients[0], optimized_coefficients[1], optimized_coefficients[2], optimized_coefficients[3], optimized_coefficients[4], optimized_coefficients[5], optimized_coefficients[6]);
335 template <
typename Po
intT,
typename Po
intNT>
void
337 const std::vector<int> &inliers,
const Eigen::VectorXf &model_coefficients,
PointCloud &projected_points,
bool copy_data_fields)
340 if (model_coefficients.size () != 7)
342 PCL_ERROR (
"[pcl::SampleConsensusModelCone::projectPoints] Invalid number of model coefficients given (%zu)!\n", model_coefficients.size ());
346 projected_points.
header = input_->header;
347 projected_points.
is_dense = input_->is_dense;
349 Eigen::Vector4f apex (model_coefficients[0], model_coefficients[1], model_coefficients[2], 0);
350 Eigen::Vector4f axis_dir (model_coefficients[3], model_coefficients[4], model_coefficients[5], 0);
351 float opening_angle = model_coefficients[6];
353 float apexdotdir = apex.dot (axis_dir);
354 float dirdotdir = 1.0f / axis_dir.dot (axis_dir);
357 if (copy_data_fields)
360 projected_points.
points.resize (input_->points.size ());
361 projected_points.
width = input_->width;
362 projected_points.
height = input_->height;
366 for (
size_t i = 0; i < projected_points.
points.size (); ++i)
371 for (
size_t i = 0; i < inliers.size (); ++i)
373 Eigen::Vector4f pt (input_->points[inliers[i]].x,
374 input_->points[inliers[i]].y,
375 input_->points[inliers[i]].z,
378 float k = (pt.dot (axis_dir) - apexdotdir) * dirdotdir;
381 pp = apex + k * axis_dir;
383 Eigen::Vector4f dir = pt - pp;
387 Eigen::Vector4f height = apex - pp;
388 float actual_cone_radius = tanf (opening_angle) * height.norm ();
391 pp += dir * actual_cone_radius;
397 projected_points.
points.resize (inliers.size ());
398 projected_points.
width =
static_cast<uint32_t
> (inliers.size ());
399 projected_points.
height = 1;
403 for (
size_t i = 0; i < inliers.size (); ++i)
408 for (
size_t i = 0; i < inliers.size (); ++i)
413 float k = (pt.dot (axis_dir) - apexdotdir) * dirdotdir;
415 pp = apex + k * axis_dir;
417 Eigen::Vector4f dir = pt - pp;
421 Eigen::Vector4f height = apex - pp;
422 float actual_cone_radius = tanf (opening_angle) * height.norm ();
425 pp += dir * actual_cone_radius;
431 template <
typename Po
intT,
typename Po
intNT>
bool
433 const std::set<int> &indices,
const Eigen::VectorXf &model_coefficients,
const double threshold)
436 if (model_coefficients.size () != 7)
438 PCL_ERROR (
"[pcl::SampleConsensusModelCone::doSamplesVerifyModel] Invalid number of model coefficients given (%zu)!\n", model_coefficients.size ());
442 Eigen::Vector4f apex (model_coefficients[0], model_coefficients[1], model_coefficients[2], 0);
443 Eigen::Vector4f axis_dir (model_coefficients[3], model_coefficients[4], model_coefficients[5], 0);
444 float openning_angle = model_coefficients[6];
446 float apexdotdir = apex.dot (axis_dir);
447 float dirdotdir = 1.0f / axis_dir.dot (axis_dir);
450 for (std::set<int>::const_iterator it = indices.begin (); it != indices.end (); ++it)
452 Eigen::Vector4f pt (input_->points[*it].x, input_->points[*it].y, input_->points[*it].z, 0);
455 float k = (pt.dot (axis_dir) - apexdotdir) * dirdotdir;
456 Eigen::Vector4f pt_proj = apex + k * axis_dir;
457 Eigen::Vector4f dir = pt - pt_proj;
461 Eigen::Vector4f height = apex - pt_proj;
462 double actual_cone_radius = tan (openning_angle) * height.norm ();
466 if (fabs (static_cast<double>(pointToAxisDistance (pt, model_coefficients) - actual_cone_radius)) > threshold)
474 template <
typename Po
intT,
typename Po
intNT>
double
476 const Eigen::Vector4f &pt,
const Eigen::VectorXf &model_coefficients)
478 Eigen::Vector4f apex (model_coefficients[0], model_coefficients[1], model_coefficients[2], 0);
479 Eigen::Vector4f axis_dir (model_coefficients[3], model_coefficients[4], model_coefficients[5], 0);
484 template <
typename Po
intT,
typename Po
intNT>
bool
488 if (model_coefficients.size () != 7)
490 PCL_ERROR (
"[pcl::SampleConsensusModelCone::isModelValid] Invalid number of model coefficients given (%zu)!\n", model_coefficients.size ());
495 if (eps_angle_ > 0.0)
498 Eigen::Vector4f coeff;
499 coeff[0] = model_coefficients[3];
500 coeff[1] = model_coefficients[4];
501 coeff[2] = model_coefficients[5];
504 Eigen::Vector4f axis (axis_[0], axis_[1], axis_[2], 0);
505 double angle_diff = fabs (
getAngle3D (axis, coeff));
506 angle_diff = (std::min) (angle_diff, M_PI - angle_diff);
508 if (angle_diff > eps_angle_)
512 if (model_coefficients[6] != -std::numeric_limits<double>::max() && model_coefficients[6] < min_angle_)
514 if (model_coefficients[6] != std::numeric_limits<double>::max() && model_coefficients[6] > max_angle_)
520 #define PCL_INSTANTIATE_SampleConsensusModelCone(PointT, PointNT) template class PCL_EXPORTS pcl::SampleConsensusModelCone<PointT, PointNT>;
522 #endif // PCL_SAMPLE_CONSENSUS_IMPL_SAC_MODEL_CONE_H_