41 #ifndef PCL_FEATURES_IMPL_ESF_H_
42 #define PCL_FEATURES_IMPL_ESF_H_
44 #include <pcl/features/esf.h>
46 #include <pcl/common/transforms.h>
50 template <
typename Po
intInT,
typename Po
intOutT>
void
54 const int binsize = 64;
55 unsigned int sample_size = 20000;
57 srand (
static_cast<unsigned int> (time (
nullptr)));
58 const auto maxindex = pc.
size ();
60 std::vector<float> d2v, d1v, d3v, wt_d3;
61 std::vector<int> wt_d2;
62 d1v.reserve (sample_size);
63 d2v.reserve (sample_size * 3);
64 d3v.reserve (sample_size);
65 wt_d2.reserve (sample_size * 3);
66 wt_d3.reserve (sample_size);
68 float h_in[binsize] = {0};
69 float h_out[binsize] = {0};
70 float h_mix[binsize] = {0};
71 float h_mix_ratio[binsize] = {0};
73 float h_a3_in[binsize] = {0};
74 float h_a3_out[binsize] = {0};
75 float h_a3_mix[binsize] = {0};
77 float h_d3_in[binsize] = {0};
78 float h_d3_out[binsize] = {0};
79 float h_d3_mix[binsize] = {0};
82 float pih =
static_cast<float>(
M_PI) / 2.0f;
86 int pcnt1,pcnt2,pcnt3;
87 for (std::size_t nn_idx = 0; nn_idx < sample_size; ++nn_idx)
90 int index1 = rand()%maxindex;
91 int index2 = rand()%maxindex;
92 int index3 = rand()%maxindex;
94 if (index1==index2 || index1 == index3 || index2 == index3)
100 Eigen::Vector4f p1 = pc[index1].getVector4fMap ();
101 Eigen::Vector4f p2 = pc[index2].getVector4fMap ();
102 Eigen::Vector4f p3 = pc[index3].getVector4fMap ();
105 Eigen::Vector4f v21 (p2 - p1);
106 Eigen::Vector4f v31 (p3 - p1);
107 Eigen::Vector4f v23 (p2 - p3);
108 a = v21.norm (); b = v31.norm (); c = v23.norm (); s = (a+b+c) * 0.5f;
109 if (s * (s-a) * (s-b) * (s-c) <= 0.001f)
120 th1 =
static_cast<int> (
pcl_round (std::acos (std::abs (v21.dot (v31))) / pih * (binsize-1)));
121 th2 =
static_cast<int> (
pcl_round (std::acos (std::abs (v23.dot (v31))) / pih * (binsize-1)));
122 th3 =
static_cast<int> (
pcl_round (std::acos (std::abs (v23.dot (v21))) / pih * (binsize-1)));
123 if (th1 < 0 || th1 >= binsize)
128 if (th2 < 0 || th2 >= binsize)
133 if (th3 < 0 || th3 >= binsize)
148 const int xs = p1[0] < 0.0?
static_cast<int>(std::floor(p1[0])+GRIDSIZE_H):
static_cast<int>(std::ceil(p1[0])+GRIDSIZE_H-1);
149 const int ys = p1[1] < 0.0?
static_cast<int>(std::floor(p1[1])+GRIDSIZE_H):
static_cast<int>(std::ceil(p1[1])+GRIDSIZE_H-1);
150 const int zs = p1[2] < 0.0?
static_cast<int>(std::floor(p1[2])+GRIDSIZE_H):
static_cast<int>(std::ceil(p1[2])+GRIDSIZE_H-1);
151 const int xt = p2[0] < 0.0?
static_cast<int>(std::floor(p2[0])+GRIDSIZE_H):
static_cast<int>(std::ceil(p2[0])+GRIDSIZE_H-1);
152 const int yt = p2[1] < 0.0?
static_cast<int>(std::floor(p2[1])+GRIDSIZE_H):
static_cast<int>(std::ceil(p2[1])+GRIDSIZE_H-1);
153 const int zt = p2[2] < 0.0?
static_cast<int>(std::floor(p2[2])+GRIDSIZE_H):
static_cast<int>(std::ceil(p2[2])+GRIDSIZE_H-1);
154 wt_d2.push_back (this->lci (xs, ys, zs, xt, yt, zt, ratio, vxlcnt, pcnt1));
155 if (wt_d2.back () == 2)
156 h_mix_ratio[
static_cast<int> (
pcl_round (ratio * (binsize-1)))]++;
157 vxlcnt_sum += vxlcnt;
162 const int xs = p1[0] < 0.0?
static_cast<int>(std::floor(p1[0])+GRIDSIZE_H):
static_cast<int>(std::ceil(p1[0])+GRIDSIZE_H-1);
163 const int ys = p1[1] < 0.0?
static_cast<int>(std::floor(p1[1])+GRIDSIZE_H):
static_cast<int>(std::ceil(p1[1])+GRIDSIZE_H-1);
164 const int zs = p1[2] < 0.0?
static_cast<int>(std::floor(p1[2])+GRIDSIZE_H):
static_cast<int>(std::ceil(p1[2])+GRIDSIZE_H-1);
165 const int xt = p3[0] < 0.0?
static_cast<int>(std::floor(p3[0])+GRIDSIZE_H):
static_cast<int>(std::ceil(p3[0])+GRIDSIZE_H-1);
166 const int yt = p3[1] < 0.0?
static_cast<int>(std::floor(p3[1])+GRIDSIZE_H):
static_cast<int>(std::ceil(p3[1])+GRIDSIZE_H-1);
167 const int zt = p3[2] < 0.0?
static_cast<int>(std::floor(p3[2])+GRIDSIZE_H):
static_cast<int>(std::ceil(p3[2])+GRIDSIZE_H-1);
168 wt_d2.push_back (this->lci (xs, ys, zs, xt, yt, zt, ratio, vxlcnt, pcnt2));
169 if (wt_d2.back () == 2)
170 h_mix_ratio[
static_cast<int>(
pcl_round (ratio * (binsize-1)))]++;
171 vxlcnt_sum += vxlcnt;
176 const int xs = p2[0] < 0.0?
static_cast<int>(std::floor(p2[0])+GRIDSIZE_H):
static_cast<int>(std::ceil(p2[0])+GRIDSIZE_H-1);
177 const int ys = p2[1] < 0.0?
static_cast<int>(std::floor(p2[1])+GRIDSIZE_H):
static_cast<int>(std::ceil(p2[1])+GRIDSIZE_H-1);
178 const int zs = p2[2] < 0.0?
static_cast<int>(std::floor(p2[2])+GRIDSIZE_H):
static_cast<int>(std::ceil(p2[2])+GRIDSIZE_H-1);
179 const int xt = p3[0] < 0.0?
static_cast<int>(std::floor(p3[0])+GRIDSIZE_H):
static_cast<int>(std::ceil(p3[0])+GRIDSIZE_H-1);
180 const int yt = p3[1] < 0.0?
static_cast<int>(std::floor(p3[1])+GRIDSIZE_H):
static_cast<int>(std::ceil(p3[1])+GRIDSIZE_H-1);
181 const int zt = p3[2] < 0.0?
static_cast<int>(std::floor(p3[2])+GRIDSIZE_H):
static_cast<int>(std::ceil(p3[2])+GRIDSIZE_H-1);
182 wt_d2.push_back (this->lci (xs,ys,zs,xt,yt,zt,ratio,vxlcnt,pcnt3));
183 if (wt_d2.back () == 2)
184 h_mix_ratio[
static_cast<int>(
pcl_round(ratio * (binsize-1)))]++;
185 vxlcnt_sum += vxlcnt;
190 d3v.push_back (std::sqrt (std::sqrt (s * (s-a) * (s-b) * (s-c))));
191 if (vxlcnt_sum <= 21)
194 h_a3_out[th1] +=
static_cast<float> (pcnt3) / 32.0f;
195 h_a3_out[th2] +=
static_cast<float> (pcnt1) / 32.0f;
196 h_a3_out[th3] +=
static_cast<float> (pcnt2) / 32.0f;
199 if (p_cnt - vxlcnt_sum < 4)
201 h_a3_in[th1] +=
static_cast<float> (pcnt3) / 32.0f;
202 h_a3_in[th2] +=
static_cast<float> (pcnt1) / 32.0f;
203 h_a3_in[th3] +=
static_cast<float> (pcnt2) / 32.0f;
208 h_a3_mix[th1] +=
static_cast<float> (pcnt3) / 32.0f;
209 h_a3_mix[th2] +=
static_cast<float> (pcnt1) / 32.0f;
210 h_a3_mix[th3] +=
static_cast<float> (pcnt2) / 32.0f;
211 wt_d3.push_back (
static_cast<float> (vxlcnt_sum) /
static_cast<float> (p_cnt));
218 for (std::size_t nn_idx = 0; nn_idx < sample_size; ++nn_idx)
221 if (d2v[nn_idx] > maxd2)
223 if (d2v[sample_size + nn_idx] > maxd2)
224 maxd2 = d2v[sample_size + nn_idx];
225 if (d2v[sample_size*2 +nn_idx] > maxd2)
226 maxd2 = d2v[sample_size*2 +nn_idx];
227 if (d3v[nn_idx] > maxd3)
233 for (std::size_t nn_idx = 0; nn_idx < sample_size; ++nn_idx)
235 if (wt_d3[nn_idx] >= 0.999)
237 index =
static_cast<int>(
pcl_round (d3v[nn_idx] / maxd3 * (binsize-1)));
238 if (index >= 0 && index < binsize)
243 if (wt_d3[nn_idx] <= 0.001)
245 index =
static_cast<int>(
pcl_round (d3v[nn_idx] / maxd3 * (binsize-1)));
246 if (index >= 0 && index < binsize)
251 index =
static_cast<int>(
pcl_round (d3v[nn_idx] / maxd3 * (binsize-1)));
252 if (index >= 0 && index < binsize)
258 for (std::size_t nn_idx = 0; nn_idx < d2v.size(); ++nn_idx )
260 if (wt_d2[nn_idx] == 0)
261 h_in[
static_cast<int>(
pcl_round (d2v[nn_idx] / maxd2 * (binsize-1)))]++ ;
262 if (wt_d2[nn_idx] == 1)
263 h_out[
static_cast<int>(
pcl_round (d2v[nn_idx] / maxd2 * (binsize-1)))]++;
264 if (wt_d2[nn_idx] == 2)
265 h_mix[
static_cast<int>(
pcl_round (d2v[nn_idx] / maxd2 * (binsize-1)))]++ ;
269 float weights[10] = {0.5f, 0.5f, 0.5f, 0.5f, 0.5f, 1.0f, 1.0f, 2.0f, 2.0f, 2.0f};
271 hist.reserve (binsize * 10);
272 for (
const float &i : h_a3_in)
273 hist.push_back (i * weights[0]);
274 for (
const float &i : h_a3_out)
275 hist.push_back (i * weights[1]);
276 for (
const float &i : h_a3_mix)
277 hist.push_back (i * weights[2]);
279 for (
const float &i : h_d3_in)
280 hist.push_back (i * weights[3]);
281 for (
const float &i : h_d3_out)
282 hist.push_back (i * weights[4]);
283 for (
const float &i : h_d3_mix)
284 hist.push_back (i * weights[5]);
286 for (
const float &i : h_in)
287 hist.push_back (i*0.5f * weights[6]);
288 for (
const float &i : h_out)
289 hist.push_back (i * weights[7]);
290 for (
const float &i : h_mix)
291 hist.push_back (i * weights[8]);
292 for (
const float &i : h_mix_ratio)
293 hist.push_back (i*0.5f * weights[9]);
296 for (
const float &i : hist)
299 for (
float &i : hist)
304 template <
typename Po
intInT,
typename Po
intOutT>
int
306 const int x1,
const int y1,
const int z1,
307 const int x2,
const int y2,
const int z2,
308 float &ratio,
int &incnt,
int &pointcount)
316 int x_inc, y_inc, z_inc;
324 int l = std::abs (dx);
329 int m = std::abs (dy);
334 int n = std::abs (dz);
338 if ((l >= m) & (l >= n))
342 for (
int i = 1; i<l; i++)
345 voxel_in +=
static_cast<int>(lut_[act_voxel[0]][act_voxel[1]][act_voxel[2]] == 1);
348 act_voxel[1] += y_inc;
353 act_voxel[2] += z_inc;
358 act_voxel[0] += x_inc;
361 else if ((m >= l) & (m >= n))
365 for (
int i=1; i<m; i++)
368 voxel_in +=
static_cast<int>(lut_[act_voxel[0]][act_voxel[1]][act_voxel[2]] == 1);
371 act_voxel[0] += x_inc;
376 act_voxel[2] += z_inc;
381 act_voxel[1] += y_inc;
388 for (
int i=1; i<n; i++)
391 voxel_in +=
static_cast<int>(lut_[act_voxel[0]][act_voxel[1]][act_voxel[2]] == 1);
394 act_voxel[1] += y_inc;
399 act_voxel[0] += x_inc;
404 act_voxel[2] += z_inc;
408 voxel_in +=
static_cast<int>(lut_[act_voxel[0]][act_voxel[1]][act_voxel[2]] == 1);
410 pointcount = voxelcount;
412 if (voxel_in >= voxelcount-1)
418 ratio =
static_cast<float>(voxel_in) /
static_cast<float>(voxelcount);
423 template <
typename Po
intInT,
typename Po
intOutT>
void
426 for (
const auto& point: cluster)
428 int xx = point.x<0.0?
static_cast<int>(std::floor(point.x)+GRIDSIZE_H) :
static_cast<int>(std::ceil(point.x)+GRIDSIZE_H-1);
429 int yy = point.y<0.0?
static_cast<int>(std::floor(point.y)+GRIDSIZE_H) :
static_cast<int>(std::ceil(point.y)+GRIDSIZE_H-1);
430 int zz = point.z<0.0?
static_cast<int>(std::floor(point.z)+GRIDSIZE_H) :
static_cast<int>(std::ceil(point.z)+GRIDSIZE_H-1);
432 for (
int x = -1; x < 2; x++)
433 for (
int y = -1; y < 2; y++)
434 for (
int z = -1; z < 2; z++)
440 if (yi >= GRIDSIZE || xi >= GRIDSIZE || zi>=GRIDSIZE || yi < 0 || xi < 0 || zi < 0)
445 this->lut_[xi][yi][zi] = 1;
451 template <
typename Po
intInT,
typename Po
intOutT>
void
454 for (
const auto& point: cluster)
456 int xx = point.x<0.0?
static_cast<int>(std::floor(point.x)+GRIDSIZE_H) :
static_cast<int>(std::ceil(point.x)+GRIDSIZE_H-1);
457 int yy = point.y<0.0?
static_cast<int>(std::floor(point.y)+GRIDSIZE_H) :
static_cast<int>(std::ceil(point.y)+GRIDSIZE_H-1);
458 int zz = point.z<0.0?
static_cast<int>(std::floor(point.z)+GRIDSIZE_H) :
static_cast<int>(std::ceil(point.z)+GRIDSIZE_H-1);
460 for (
int x = -1; x < 2; x++)
461 for (
int y = -1; y < 2; y++)
462 for (
int z = -1; z < 2; z++)
468 if (yi >= GRIDSIZE || xi >= GRIDSIZE || zi>=GRIDSIZE || yi < 0 || xi < 0 || zi < 0)
473 this->lut_[xi][yi][zi] = 0;
479 template <
typename Po
intInT,
typename Po
intOutT>
void
486 float max_distance = 0;
489 for (
const auto& point: local_cloud_)
492 if (d > max_distance)
496 float scale_factor = 1.0f / max_distance * scalefactor;
498 Eigen::Affine3f matrix = Eigen::Affine3f::Identity();
499 matrix.scale (scale_factor);
504 template<
typename Po
intInT,
typename Po
intOutT>
void
509 output.width = output.height = 0;
510 output.points.clear ();
514 output.header = input_->header;
520 output.width = output.height = 1;
521 output.is_dense = input_->is_dense;
522 output.points.resize (1);
525 computeFeature (output);
532 template <
typename Po
intInT,
typename Po
intOutT>
void
535 Eigen::Vector4f xyz_centroid;
536 std::vector<float> hist;
537 scale_points_unit_sphere (*surface_,
static_cast<float>(GRIDSIZE_H), xyz_centroid);
538 this->voxelize9 (local_cloud_);
539 this->computeESF (local_cloud_, hist);
540 this->cleanup9 (local_cloud_);
543 output.points.resize (1);
547 for (std::size_t d = 0; d < hist.size (); ++d)
548 output[0].histogram[d] = hist[d];
551 #define PCL_INSTANTIATE_ESFEstimation(T,OutT) template class PCL_EXPORTS pcl::ESFEstimation<T,OutT>;
553 #endif // PCL_FEATURES_IMPL_ESF_H_