41 #ifndef PCL_SAMPLE_CONSENSUS_IMPL_SAC_MODEL_SPHERE_H_
42 #define PCL_SAMPLE_CONSENSUS_IMPL_SAC_MODEL_SPHERE_H_
44 #include <pcl/sample_consensus/sac_model_sphere.h>
47 template <
typename Po
intT>
bool
50 if (samples.size () != sample_size_)
52 PCL_ERROR (
"[pcl::SampleConsensusModelSphere::isSampleGood] Wrong number of samples (is %lu, should be %lu)!\n", samples.size (), sample_size_);
59 template <
typename Po
intT>
bool
61 const Indices &samples, Eigen::VectorXf &model_coefficients)
const
64 if (samples.size () != sample_size_)
66 PCL_ERROR (
"[pcl::SampleConsensusModelSphere::computeModelCoefficients] Invalid set of samples given (%lu)!\n", samples.size ());
72 for (
int i = 0; i < 4; i++)
74 temp (i, 0) = (*input_)[samples[i]].x;
75 temp (i, 1) = (*input_)[samples[i]].y;
76 temp (i, 2) = (*input_)[samples[i]].z;
79 const double m11 = temp.determinant ();
85 for (
int i = 0; i < 4; ++i)
87 temp (i, 0) = ((*input_)[samples[i]].x) * ((*input_)[samples[i]].x) +
88 ((*input_)[samples[i]].y) * ((*input_)[samples[i]].y) +
89 ((*input_)[samples[i]].z) * ((*input_)[samples[i]].z);
91 const double m12 = temp.determinant ();
93 for (
int i = 0; i < 4; ++i)
95 temp (i, 1) = temp (i, 0);
96 temp (i, 0) = (*input_)[samples[i]].x;
98 const double m13 = temp.determinant ();
100 for (
int i = 0; i < 4; ++i)
102 temp (i, 2) = temp (i, 1);
103 temp (i, 1) = (*input_)[samples[i]].y;
105 const double m14 = temp.determinant ();
107 for (
int i = 0; i < 4; ++i)
109 temp (i, 0) = temp (i, 2);
110 temp (i, 1) = (*input_)[samples[i]].x;
111 temp (i, 2) = (*input_)[samples[i]].y;
112 temp (i, 3) = (*input_)[samples[i]].z;
114 const double m15 = temp.determinant ();
117 model_coefficients.resize (model_size_);
118 model_coefficients[0] = 0.5f * m12 / m11;
119 model_coefficients[1] = 0.5f * m13 / m11;
120 model_coefficients[2] = 0.5f * m14 / m11;
122 model_coefficients[3] = std::sqrt (model_coefficients[0] * model_coefficients[0] +
123 model_coefficients[1] * model_coefficients[1] +
124 model_coefficients[2] * model_coefficients[2] - m15 / m11);
126 PCL_DEBUG (
"[pcl::SampleConsensusModelSphere::computeModelCoefficients] Model is (%g,%g,%g,%g)\n",
127 model_coefficients[0], model_coefficients[1], model_coefficients[2], model_coefficients[3]);
131 #define AT(POS) ((*input_)[(*indices_)[(POS)]])
137 const __m256 tmp1 = _mm256_sub_ps (_mm256_set_ps (AT(i ).x, AT(i+1).x, AT(i+2).x, AT(i+3).x, AT(i+4).x, AT(i+5).x, AT(i+6).x, AT(i+7).x), a_vec);
138 const __m256 tmp2 = _mm256_sub_ps (_mm256_set_ps (AT(i ).y, AT(i+1).y, AT(i+2).y, AT(i+3).y, AT(i+4).y, AT(i+5).y, AT(i+6).y, AT(i+7).y), b_vec);
139 const __m256 tmp3 = _mm256_sub_ps (_mm256_set_ps (AT(i ).z, AT(i+1).z, AT(i+2).z, AT(i+3).z, AT(i+4).z, AT(i+5).z, AT(i+6).z, AT(i+7).z), c_vec);
140 return _mm256_add_ps (_mm256_add_ps (_mm256_mul_ps (tmp1, tmp1), _mm256_mul_ps (tmp2, tmp2)), _mm256_mul_ps(tmp3, tmp3));
148 const __m128 tmp1 = _mm_sub_ps (_mm_set_ps (AT(i ).x, AT(i+1).x, AT(i+2).x, AT(i+3).x), a_vec);
149 const __m128 tmp2 = _mm_sub_ps (_mm_set_ps (AT(i ).y, AT(i+1).y, AT(i+2).y, AT(i+3).y), b_vec);
150 const __m128 tmp3 = _mm_sub_ps (_mm_set_ps (AT(i ).z, AT(i+1).z, AT(i+2).z, AT(i+3).z), c_vec);
151 return _mm_add_ps (_mm_add_ps (_mm_mul_ps (tmp1, tmp1), _mm_mul_ps (tmp2, tmp2)), _mm_mul_ps(tmp3, tmp3));
158 template <
typename Po
intT>
void
160 const Eigen::VectorXf &model_coefficients, std::vector<double> &distances)
const
163 if (!isModelValid (model_coefficients))
168 distances.resize (indices_->size ());
170 const Eigen::Vector3f center (model_coefficients[0], model_coefficients[1], model_coefficients[2]);
172 for (std::size_t i = 0; i < indices_->size (); ++i)
176 distances[i] = std::abs (((*input_)[(*indices_)[i]].getVector3fMap () - center).norm () - model_coefficients[3]);
181 template <
typename Po
intT>
void
183 const Eigen::VectorXf &model_coefficients,
const double threshold,
Indices &inliers)
186 if (!isModelValid (model_coefficients))
193 error_sqr_dists_.clear ();
194 inliers.reserve (indices_->size ());
195 error_sqr_dists_.reserve (indices_->size ());
197 const float sqr_inner_radius = (model_coefficients[3] <= threshold ? 0.0f : (model_coefficients[3] - threshold) * (model_coefficients[3] - threshold));
198 const float sqr_outer_radius = (model_coefficients[3] + threshold) * (model_coefficients[3] + threshold);
199 const Eigen::Vector3f center (model_coefficients[0], model_coefficients[1], model_coefficients[2]);
201 for (std::size_t i = 0; i < indices_->size (); ++i)
205 const float sqr_dist = ((*input_)[(*indices_)[i]].getVector3fMap () - center).squaredNorm ();
206 if ((sqr_dist <= sqr_outer_radius) && (sqr_dist >= sqr_inner_radius))
209 inliers.push_back ((*indices_)[i]);
211 error_sqr_dists_.push_back (
static_cast<double> (std::abs (std::sqrt (sqr_dist) - model_coefficients[3])));
217 template <
typename Po
intT> std::size_t
219 const Eigen::VectorXf &model_coefficients,
const double threshold)
const
222 if (!isModelValid (model_coefficients))
225 #if defined (__AVX__) && defined (__AVX2__)
226 return countWithinDistanceAVX (model_coefficients, threshold);
227 #elif defined (__SSE__) && defined (__SSE2__) && defined (__SSE4_1__)
228 return countWithinDistanceSSE (model_coefficients, threshold);
230 return countWithinDistanceStandard (model_coefficients, threshold);
235 template <
typename Po
intT> std::size_t
237 const Eigen::VectorXf &model_coefficients,
const double threshold, std::size_t i)
const
239 std::size_t nr_p = 0;
240 const float sqr_inner_radius = (model_coefficients[3] <= threshold ? 0.0f : (model_coefficients[3] - threshold) * (model_coefficients[3] - threshold));
241 const float sqr_outer_radius = (model_coefficients[3] + threshold) * (model_coefficients[3] + threshold);
242 const Eigen::Vector3f center (model_coefficients[0], model_coefficients[1], model_coefficients[2]);
244 for (; i < indices_->size (); ++i)
248 const float sqr_dist = ((*input_)[(*indices_)[i]].getVector3fMap () - center).squaredNorm ();
249 if ((sqr_dist <= sqr_outer_radius) && (sqr_dist >= sqr_inner_radius))
256 #if defined (__SSE__) && defined (__SSE2__) && defined (__SSE4_1__)
257 template <
typename Po
intT> std::size_t
259 const Eigen::VectorXf &model_coefficients,
const double threshold, std::size_t i)
const
261 std::size_t nr_p = 0;
262 const __m128 a_vec = _mm_set1_ps (model_coefficients[0]);
263 const __m128 b_vec = _mm_set1_ps (model_coefficients[1]);
264 const __m128 c_vec = _mm_set1_ps (model_coefficients[2]);
266 const __m128 sqr_inner_radius = _mm_set1_ps ((model_coefficients[3] <= threshold ? 0.0 : (model_coefficients[3]-threshold)*(model_coefficients[3]-threshold)));
267 const __m128 sqr_outer_radius = _mm_set1_ps ((model_coefficients[3]+threshold)*(model_coefficients[3]+threshold));
268 __m128i res = _mm_set1_epi32(0);
269 for (; (i + 4) <= indices_->size (); i += 4)
271 const __m128 sqr_dist = sqr_dist4 (i, a_vec, b_vec, c_vec);
272 const __m128 mask = _mm_and_ps (_mm_cmplt_ps (sqr_inner_radius, sqr_dist), _mm_cmplt_ps (sqr_dist, sqr_outer_radius));
273 res = _mm_add_epi32 (res, _mm_and_si128 (_mm_set1_epi32 (1), _mm_castps_si128 (mask)));
280 nr_p += _mm_extract_epi32 (res, 0);
281 nr_p += _mm_extract_epi32 (res, 1);
282 nr_p += _mm_extract_epi32 (res, 2);
283 nr_p += _mm_extract_epi32 (res, 3);
286 nr_p += countWithinDistanceStandard (model_coefficients, threshold, i);
292 #if defined (__AVX__) && defined (__AVX2__)
293 template <
typename Po
intT> std::size_t
295 const Eigen::VectorXf &model_coefficients,
const double threshold, std::size_t i)
const
297 std::size_t nr_p = 0;
298 const __m256 a_vec = _mm256_set1_ps (model_coefficients[0]);
299 const __m256 b_vec = _mm256_set1_ps (model_coefficients[1]);
300 const __m256 c_vec = _mm256_set1_ps (model_coefficients[2]);
302 const __m256 sqr_inner_radius = _mm256_set1_ps ((model_coefficients[3] <= threshold ? 0.0 : (model_coefficients[3]-threshold)*(model_coefficients[3]-threshold)));
303 const __m256 sqr_outer_radius = _mm256_set1_ps ((model_coefficients[3]+threshold)*(model_coefficients[3]+threshold));
304 __m256i res = _mm256_set1_epi32(0);
305 for (; (i + 8) <= indices_->size (); i += 8)
307 const __m256 sqr_dist = sqr_dist8 (i, a_vec, b_vec, c_vec);
308 const __m256 mask = _mm256_and_ps (_mm256_cmp_ps (sqr_inner_radius, sqr_dist, _CMP_LT_OQ), _mm256_cmp_ps (sqr_dist, sqr_outer_radius, _CMP_LT_OQ));
309 res = _mm256_add_epi32 (res, _mm256_and_si256 (_mm256_set1_epi32 (1), _mm256_castps_si256 (mask)));
320 nr_p += _mm256_extract_epi32 (res, 0);
321 nr_p += _mm256_extract_epi32 (res, 1);
322 nr_p += _mm256_extract_epi32 (res, 2);
323 nr_p += _mm256_extract_epi32 (res, 3);
324 nr_p += _mm256_extract_epi32 (res, 4);
325 nr_p += _mm256_extract_epi32 (res, 5);
326 nr_p += _mm256_extract_epi32 (res, 6);
327 nr_p += _mm256_extract_epi32 (res, 7);
330 nr_p += countWithinDistanceStandard (model_coefficients, threshold, i);
336 template <
typename Po
intT>
void
338 const Indices &inliers,
const Eigen::VectorXf &model_coefficients, Eigen::VectorXf &optimized_coefficients)
const
340 optimized_coefficients = model_coefficients;
343 if (!isModelValid (model_coefficients))
345 PCL_ERROR (
"[pcl::SampleConsensusModelSphere::optimizeModelCoefficients] Given model is invalid!\n");
350 if (inliers.size () <= sample_size_)
352 PCL_ERROR (
"[pcl::SampleConsensusModelSphere::optimizeModelCoefficients] Not enough inliers to refine/optimize the model's coefficients (%lu)! Returning the same coefficients.\n", inliers.size ());
356 Eigen::ArrayXf pts_x(inliers.size());
357 Eigen::ArrayXf pts_y(inliers.size());
358 Eigen::ArrayXf pts_z(inliers.size());
360 for(
const auto& index : inliers) {
361 pts_x[pos] = (*input_)[index].x;
362 pts_y[pos] = (*input_)[index].y;
363 pts_z[pos] = (*input_)[index].z;
368 PCL_DEBUG (
"[pcl::SampleConsensusModelSphere::optimizeModelCoefficients] Initial solution: %g %g %g %g \nFinal solution: %g %g %g %g\n",
369 model_coefficients[0], model_coefficients[1], model_coefficients[2], model_coefficients[3], optimized_coefficients[0], optimized_coefficients[1], optimized_coefficients[2], optimized_coefficients[3]);
373 template <
typename Po
intT>
void
375 const Indices &inliers,
const Eigen::VectorXf &model_coefficients,
PointCloud &projected_points,
bool copy_data_fields)
const
378 if (!isModelValid (model_coefficients))
380 PCL_ERROR (
"[pcl::SampleConsensusModelSphere::projectPoints] Given model is invalid!\n");
384 projected_points.
header = input_->header;
385 projected_points.
is_dense = input_->is_dense;
388 const Eigen::Vector3d C (model_coefficients[0], model_coefficients[1], model_coefficients[2]);
390 const double r = model_coefficients[3];
393 if (copy_data_fields)
396 projected_points.
resize (input_->size ());
397 projected_points.
width = input_->width;
398 projected_points.
height = input_->height;
400 using FieldList =
typename pcl::traits::fieldList<PointT>::type;
402 for (std::size_t i = 0; i < projected_points.
points.size (); ++i)
407 for (
const auto& inlier : inliers)
411 const Eigen::Vector3d P (input_->points[inlier].x, input_->points[inlier].y, input_->points[inlier].z);
413 const Eigen::Vector3d direction = (P - C).normalized();
416 const Eigen::Vector3d
K = C + r * direction;
418 projected_points.
points[inlier].x =
static_cast<float> (
K[0]);
419 projected_points.
points[inlier].y =
static_cast<float> (
K[1]);
420 projected_points.
points[inlier].z =
static_cast<float> (
K[2]);
426 projected_points.
resize (inliers.size ());
427 projected_points.
width =
static_cast<uint32_t
> (inliers.size ());
428 projected_points.
height = 1;
430 using FieldList =
typename pcl::traits::fieldList<PointT>::type;
432 for (std::size_t i = 0; i < inliers.size (); ++i)
437 for (std::size_t i = 0; i < inliers.size (); ++i)
441 const Eigen::Vector3d P (input_->points[inliers[i]].x, input_->points[inliers[i]].y, input_->points[inliers[i]].z);
443 const Eigen::Vector3d direction = (P - C).normalized();
446 const Eigen::Vector3d
K = C + r * direction;
448 projected_points.
points[i].x =
static_cast<float> (
K[0]);
449 projected_points.
points[i].y =
static_cast<float> (
K[1]);
450 projected_points.
points[i].z =
static_cast<float> (
K[2]);
456 template <
typename Po
intT>
bool
458 const std::set<index_t> &indices,
const Eigen::VectorXf &model_coefficients,
const double threshold)
const
461 if (!isModelValid (model_coefficients))
463 PCL_ERROR (
"[pcl::SampleConsensusModelSphere::doSamplesVerifyModel] Given model is invalid!\n");
467 const float sqr_inner_radius = (model_coefficients[3] <= threshold ? 0.0f : (model_coefficients[3] - threshold) * (model_coefficients[3] - threshold));
468 const float sqr_outer_radius = (model_coefficients[3] + threshold) * (model_coefficients[3] + threshold);
469 const Eigen::Vector3f center (model_coefficients[0], model_coefficients[1], model_coefficients[2]);
470 for (
const auto &index : indices)
474 const float sqr_dist = ((*input_)[index].getVector3fMap () - center).squaredNorm ();
475 if ((sqr_dist > sqr_outer_radius) || (sqr_dist < sqr_inner_radius))
484 #define PCL_INSTANTIATE_SampleConsensusModelSphere(T) template class PCL_EXPORTS pcl::SampleConsensusModelSphere<T>;
PointCloud represents the base class in PCL for storing collections of 3D points.
bool is_dense
True if no points are invalid (e.g., have NaN or Inf values in any of their floating point fields).
void resize(std::size_t count)
Resizes the container to contain count elements.
std::uint32_t width
The point cloud width (if organized as an image-structure).
pcl::PCLHeader header
The point cloud header.
std::uint32_t height
The point cloud height (if organized as an image-structure).
std::vector< PointT, Eigen::aligned_allocator< PointT > > points
The point data.
SampleConsensusModelSphere defines a model for 3D sphere segmentation.
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 sphere model.
void optimizeModelCoefficients(const Indices &inliers, const Eigen::VectorXf &model_coefficients, Eigen::VectorXf &optimized_coefficients) const override
Recompute the sphere 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 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 sphere model coefficients.
std::size_t countWithinDistanceStandard(const Eigen::VectorXf &model_coefficients, const double threshold, std::size_t i=0) const
This implementation uses no SIMD instructions.
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 sphere model.
bool computeModelCoefficients(const Indices &samples, Eigen::VectorXf &model_coefficients) const override
Check whether the given index samples can form a valid sphere model, compute the model coefficients f...
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.
PCL_EXPORTS int optimizeModelCoefficientsSphere(Eigen::VectorXf &coeff, const Eigen::ArrayXf &pts_x, const Eigen::ArrayXf &pts_y, const Eigen::ArrayXf &pts_z)
IndicesAllocator<> Indices
Type used for indices in PCL.
Helper functor structure for concatenate.