41 #ifndef PCL_FEATURES_IMPL_ESF_H_
42 #define PCL_FEATURES_IMPL_ESF_H_
50 template <
typename Po
intInT,
typename Po
intOutT>
void
54 const int binsize = 64;
55 unsigned int sample_size = 20000;
56 srand (static_cast<unsigned int> (time (0)));
57 int maxindex =
static_cast<int> (pc.points.size ());
59 int index1, index2, index3;
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};
76 float h_d1[binsize] = {0};
78 float h_d3_in[binsize] = {0};
79 float h_d3_out[binsize] = {0};
80 float h_d3_mix[binsize] = {0};
83 float pih =
static_cast<float>(M_PI) / 2.0f;
87 int pcnt1,pcnt2,pcnt3;
88 for (
size_t nn_idx = 0; nn_idx < sample_size; ++nn_idx)
91 index1 = rand()%maxindex;
92 index2 = rand()%maxindex;
93 index3 = rand()%maxindex;
95 if (index1==index2 || index1 == index3 || index2 == index3)
101 Eigen::Vector4f p1 = pc.points[index1].getVector4fMap ();
102 Eigen::Vector4f p2 = pc.points[index2].getVector4fMap ();
103 Eigen::Vector4f p3 = pc.points[index3].getVector4fMap ();
106 Eigen::Vector4f v21 (p2 - p1);
107 Eigen::Vector4f v31 (p3 - p1);
108 Eigen::Vector4f v23 (p2 - p3);
109 a = v21.norm (); b = v31.norm (); c = v23.norm (); s = (a+b+c) * 0.5f;
110 if (s * (s-a) * (s-b) * (s-c) <= 0.001f)
118 th1 =
static_cast<int> (
pcl_round (acos (fabs (v21.dot (v31))) / pih * (binsize-1)));
119 th2 =
static_cast<int> (
pcl_round (acos (fabs (v23.dot (v31))) / pih * (binsize-1)));
120 th3 =
static_cast<int> (
pcl_round (acos (fabs (v23.dot (v21))) / pih * (binsize-1)));
121 if (th1 < 0 || th1 >= binsize)
126 if (th2 < 0 || th2 >= binsize)
131 if (th3 < 0 || th3 >= binsize)
150 const int xs = p1[0] < 0.0?
static_cast<int>(floor(p1[0])+
GRIDSIZE_H): static_cast<int>(ceil(p1[0])+
GRIDSIZE_H-1);
151 const int ys = p1[1] < 0.0?
static_cast<int>(floor(p1[1])+
GRIDSIZE_H): static_cast<int>(ceil(p1[1])+
GRIDSIZE_H-1);
152 const int zs = p1[2] < 0.0?
static_cast<int>(floor(p1[2])+
GRIDSIZE_H): static_cast<int>(ceil(p1[2])+
GRIDSIZE_H-1);
153 const int xt = p2[0] < 0.0?
static_cast<int>(floor(p2[0])+
GRIDSIZE_H): static_cast<int>(ceil(p2[0])+
GRIDSIZE_H-1);
154 const int yt = p2[1] < 0.0?
static_cast<int>(floor(p2[1])+
GRIDSIZE_H): static_cast<int>(ceil(p2[1])+
GRIDSIZE_H-1);
155 const int zt = p2[2] < 0.0?
static_cast<int>(floor(p2[2])+
GRIDSIZE_H): static_cast<int>(ceil(p2[2])+
GRIDSIZE_H-1);
156 wt_d2.push_back (this->lci (xs, ys, zs, xt, yt, zt, ratio, vxlcnt, pcnt1));
157 if (wt_d2.back () == 2)
158 h_mix_ratio[static_cast<int> (
pcl_round (ratio * (binsize-1)))]++;
159 vxlcnt_sum += vxlcnt;
164 const int xs = p1[0] < 0.0?
static_cast<int>(floor(p1[0])+
GRIDSIZE_H): static_cast<int>(ceil(p1[0])+
GRIDSIZE_H-1);
165 const int ys = p1[1] < 0.0?
static_cast<int>(floor(p1[1])+
GRIDSIZE_H): static_cast<int>(ceil(p1[1])+
GRIDSIZE_H-1);
166 const int zs = p1[2] < 0.0?
static_cast<int>(floor(p1[2])+
GRIDSIZE_H): static_cast<int>(ceil(p1[2])+
GRIDSIZE_H-1);
167 const int xt = p3[0] < 0.0?
static_cast<int>(floor(p3[0])+
GRIDSIZE_H): static_cast<int>(ceil(p3[0])+
GRIDSIZE_H-1);
168 const int yt = p3[1] < 0.0?
static_cast<int>(floor(p3[1])+
GRIDSIZE_H): static_cast<int>(ceil(p3[1])+
GRIDSIZE_H-1);
169 const int zt = p3[2] < 0.0?
static_cast<int>(floor(p3[2])+
GRIDSIZE_H): static_cast<int>(ceil(p3[2])+
GRIDSIZE_H-1);
170 wt_d2.push_back (this->lci (xs, ys, zs, xt, yt, zt, ratio, vxlcnt, pcnt2));
171 if (wt_d2.back () == 2)
172 h_mix_ratio[static_cast<int>(
pcl_round (ratio * (binsize-1)))]++;
173 vxlcnt_sum += vxlcnt;
178 const int xs = p2[0] < 0.0?
static_cast<int>(floor(p2[0])+
GRIDSIZE_H): static_cast<int>(ceil(p2[0])+
GRIDSIZE_H-1);
179 const int ys = p2[1] < 0.0?
static_cast<int>(floor(p2[1])+
GRIDSIZE_H): static_cast<int>(ceil(p2[1])+
GRIDSIZE_H-1);
180 const int zs = p2[2] < 0.0?
static_cast<int>(floor(p2[2])+
GRIDSIZE_H): static_cast<int>(ceil(p2[2])+
GRIDSIZE_H-1);
181 const int xt = p3[0] < 0.0?
static_cast<int>(floor(p3[0])+
GRIDSIZE_H): static_cast<int>(ceil(p3[0])+
GRIDSIZE_H-1);
182 const int yt = p3[1] < 0.0?
static_cast<int>(floor(p3[1])+
GRIDSIZE_H): static_cast<int>(ceil(p3[1])+
GRIDSIZE_H-1);
183 const int zt = p3[2] < 0.0?
static_cast<int>(floor(p3[2])+
GRIDSIZE_H): static_cast<int>(ceil(p3[2])+
GRIDSIZE_H-1);
184 wt_d2.push_back (this->lci (xs,ys,zs,xt,yt,zt,ratio,vxlcnt,pcnt3));
185 if (wt_d2.back () == 2)
186 h_mix_ratio[static_cast<int>(
pcl_round(ratio * (binsize-1)))]++;
187 vxlcnt_sum += vxlcnt;
192 d3v.push_back (
sqrt (
sqrt (s * (s-a) * (s-b) * (s-c))));
193 if (vxlcnt_sum <= 21)
196 h_a3_out[th1] +=
static_cast<float> (pcnt3) / 32.0f;
197 h_a3_out[th2] +=
static_cast<float> (pcnt1) / 32.0f;
198 h_a3_out[th3] +=
static_cast<float> (pcnt2) / 32.0f;
201 if (p_cnt - vxlcnt_sum < 4)
203 h_a3_in[th1] +=
static_cast<float> (pcnt3) / 32.0f;
204 h_a3_in[th2] +=
static_cast<float> (pcnt1) / 32.0f;
205 h_a3_in[th3] +=
static_cast<float> (pcnt2) / 32.0f;
210 h_a3_mix[th1] +=
static_cast<float> (pcnt3) / 32.0f;
211 h_a3_mix[th2] +=
static_cast<float> (pcnt1) / 32.0f;
212 h_a3_mix[th3] +=
static_cast<float> (pcnt2) / 32.0f;
213 wt_d3.push_back (static_cast<float> (vxlcnt_sum) / static_cast<float> (p_cnt));
221 for (
size_t nn_idx = 0; nn_idx < sample_size; ++nn_idx)
224 if (d1v[nn_idx] > maxd1)
226 if (d2v[nn_idx] > maxd2)
228 if (d2v[sample_size + nn_idx] > maxd2)
229 maxd2 = d2v[sample_size + nn_idx];
230 if (d2v[sample_size*2 +nn_idx] > maxd2)
231 maxd2 = d2v[sample_size*2 +nn_idx];
232 if (d3v[nn_idx] > maxd3)
238 for (
size_t nn_idx = 0; nn_idx < sample_size; ++nn_idx)
240 h_d1[
static_cast<int>(
pcl_round (d1v[nn_idx] / maxd1 * (binsize-1)))]++ ;
242 if (wt_d3[nn_idx] >= 0.999)
244 index =
static_cast<int>(
pcl_round (d3v[nn_idx] / maxd3 * (binsize-1)));
245 if (index >= 0 && index < binsize)
250 if (wt_d3[nn_idx] <= 0.001)
252 index =
static_cast<int>(
pcl_round (d3v[nn_idx] / maxd3 * (binsize-1)));
253 if (index >= 0 && index < binsize)
258 index =
static_cast<int>(
pcl_round (d3v[nn_idx] / maxd3 * (binsize-1)));
259 if (index >= 0 && index < binsize)
265 for (
size_t nn_idx = 0; nn_idx < d2v.size(); ++nn_idx )
267 if (wt_d2[nn_idx] == 0)
268 h_in[
static_cast<int>(
pcl_round (d2v[nn_idx] / maxd2 * (binsize-1)))]++ ;
269 if (wt_d2[nn_idx] == 1)
270 h_out[
static_cast<int>(
pcl_round (d2v[nn_idx] / maxd2 * (binsize-1)))]++;
271 if (wt_d2[nn_idx] == 2)
272 h_mix[
static_cast<int>(
pcl_round (d2v[nn_idx] / maxd2 * (binsize-1)))]++ ;
276 float weights[10] = {0.5f, 0.5f, 0.5f, 0.5f, 0.5f, 1.0f, 1.0f, 2.0f, 2.0f, 2.0f};
278 hist.reserve (binsize * 10);
279 for (
int i = 0; i < binsize; i++)
280 hist.push_back (h_a3_in[i] * weights[0]);
281 for (
int i = 0; i < binsize; i++)
282 hist.push_back (h_a3_out[i] * weights[1]);
283 for (
int i = 0; i < binsize; i++)
284 hist.push_back (h_a3_mix[i] * weights[2]);
286 for (
int i = 0; i < binsize; i++)
287 hist.push_back (h_d3_in[i] * weights[3]);
288 for (
int i = 0; i < binsize; i++)
289 hist.push_back (h_d3_out[i] * weights[4]);
290 for (
int i = 0; i < binsize; i++)
291 hist.push_back (h_d3_mix[i] * weights[5]);
293 for (
int i = 0; i < binsize; i++)
294 hist.push_back (h_in[i]*0.5f * weights[6]);
295 for (
int i = 0; i < binsize; i++)
296 hist.push_back (h_out[i] * weights[7]);
297 for (
int i = 0; i < binsize; i++)
298 hist.push_back (h_mix[i] * weights[8]);
299 for (
int i = 0; i < binsize; i++)
300 hist.push_back (h_mix_ratio[i]*0.5f * weights[9]);
303 for (
size_t i = 0; i < hist.size (); i++)
306 for (
size_t i = 0; i < hist.size (); i++)
311 template <
typename Po
intInT,
typename Po
intOutT>
int
313 const int x1,
const int y1,
const int z1,
314 const int x2,
const int y2,
const int z2,
315 float &ratio,
int &incnt,
int &pointcount)
323 int x_inc, y_inc, z_inc;
345 if ((l >= m) & (l >= n))
349 for (
int i = 1; i<l; i++)
352 voxel_in +=
static_cast<int>(lut_[act_voxel[0]][act_voxel[1]][act_voxel[2]] == 1);
355 act_voxel[1] += y_inc;
360 act_voxel[2] += z_inc;
365 act_voxel[0] += x_inc;
368 else if ((m >= l) & (m >= n))
372 for (
int i=1; i<m; i++)
375 voxel_in +=
static_cast<int>(lut_[act_voxel[0]][act_voxel[1]][act_voxel[2]] == 1);
378 act_voxel[0] += x_inc;
383 act_voxel[2] += z_inc;
388 act_voxel[1] += y_inc;
395 for (
int i=1; i<n; i++)
398 voxel_in +=
static_cast<int>(lut_[act_voxel[0]][act_voxel[1]][act_voxel[2]] == 1);
401 act_voxel[1] += y_inc;
406 act_voxel[0] += x_inc;
411 act_voxel[2] += z_inc;
415 voxel_in +=
static_cast<int>(lut_[act_voxel[0]][act_voxel[1]][act_voxel[2]] == 1);
417 pointcount = voxelcount;
419 if (voxel_in >= voxelcount-1)
425 ratio =
static_cast<float>(voxel_in) / static_cast<float>(voxelcount);
430 template <
typename Po
intInT,
typename Po
intOutT>
void
433 int xi,yi,zi,xx,yy,zz;
434 for (
size_t i = 0; i < cluster.points.size (); ++i)
436 xx = cluster.points[i].x<0.0?
static_cast<int>(floor(cluster.points[i].x)+
GRIDSIZE_H) : static_cast<int>(ceil(cluster.points[i].x)+
GRIDSIZE_H-1);
437 yy = cluster.points[i].y<0.0?
static_cast<int>(floor(cluster.points[i].y)+
GRIDSIZE_H) : static_cast<int>(ceil(cluster.points[i].y)+
GRIDSIZE_H-1);
438 zz = cluster.points[i].z<0.0?
static_cast<int>(floor(cluster.points[i].z)+
GRIDSIZE_H) : static_cast<int>(ceil(cluster.points[i].z)+
GRIDSIZE_H-1);
440 for (
int x = -1; x < 2; x++)
441 for (
int y = -1; y < 2; y++)
442 for (
int z = -1; z < 2; z++)
453 this->lut_[xi][yi][zi] = 1;
459 template <
typename Po
intInT,
typename Po
intOutT>
void
462 int xi,yi,zi,xx,yy,zz;
463 for (
size_t i = 0; i < cluster.points.size (); ++i)
465 xx = cluster.points[i].x<0.0?
static_cast<int>(floor(cluster.points[i].x)+
GRIDSIZE_H) : static_cast<int>(ceil(cluster.points[i].x)+
GRIDSIZE_H-1);
466 yy = cluster.points[i].y<0.0?
static_cast<int>(floor(cluster.points[i].y)+
GRIDSIZE_H) : static_cast<int>(ceil(cluster.points[i].y)+
GRIDSIZE_H-1);
467 zz = cluster.points[i].z<0.0?
static_cast<int>(floor(cluster.points[i].z)+
GRIDSIZE_H) : static_cast<int>(ceil(cluster.points[i].z)+
GRIDSIZE_H-1);
469 for (
int x = -1; x < 2; x++)
470 for (
int y = -1; y < 2; y++)
471 for (
int z = -1; z < 2; z++)
482 this->lut_[xi][yi][zi] = 0;
488 template <
typename Po
intInT,
typename Po
intOutT>
void
495 float max_distance = 0, d;
498 for (
size_t i = 0; i < local_cloud_.points.size (); ++i)
501 if (d > max_distance)
505 float scale_factor = 1.0f / max_distance * scalefactor;
507 Eigen::Affine3f matrix = Eigen::Affine3f::Identity();
508 matrix.scale (scale_factor);
513 template<
typename Po
intInT,
typename Po
intOutT>
void
523 output.
header = input_->header;
534 computeFeature (output);
541 template <
typename Po
intInT,
typename Po
intOutT>
void
544 Eigen::Vector4f xyz_centroid;
545 std::vector<float> hist;
546 scale_points_unit_sphere (*surface_, static_cast<float>(
GRIDSIZE_H), xyz_centroid);
547 this->voxelize9 (local_cloud_);
548 this->computeESF (local_cloud_, hist);
549 this->cleanup9 (local_cloud_);
552 output.points.resize (1);
556 for (
size_t d = 0; d < hist.size (); ++d)
557 output.points[0].histogram[d] = hist[d];
560 #define PCL_INSTANTIATE_ESFEstimation(T,OutT) template class PCL_EXPORTS pcl::ESFEstimation<T,OutT>;
562 #endif // PCL_FEATURES_IMPL_ESF_H_