39#ifndef PCL_OCTREE_SEARCH_IMPL_H_
40#define PCL_OCTREE_SEARCH_IMPL_H_
48template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
51 const PointT& point,
Indices& point_idx_data)
54 "Invalid (NaN, Inf) point coordinates given to nearestKSearch!");
56 bool b_success =
false;
61 LeafContainerT* leaf = this->
findLeaf(key);
64 (*leaf).getPointIndices(point_idx_data);
71template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
77 return (this->
voxelSearch(search_point, point_idx_data));
80template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
86 std::vector<float>& k_sqr_distances)
90 "Invalid (NaN, Inf) point coordinates given to nearestKSearch!");
93 k_sqr_distances.clear();
98 prioPointQueueEntry point_entry;
99 std::vector<prioPointQueueEntry> point_candidates;
100 point_candidates.reserve(k);
103 key.
x = key.
y = key.
z = 0;
106 double smallest_dist = std::numeric_limits<double>::max();
109 p_q, k, this->
root_node_, key, 1, smallest_dist, point_candidates);
111 const auto result_count =
static_cast<uindex_t>(point_candidates.size());
113 k_indices.resize(result_count);
114 k_sqr_distances.resize(result_count);
116 for (
uindex_t i = 0; i < result_count; ++i) {
117 k_indices[i] = point_candidates[i].point_idx_;
118 k_sqr_distances[i] = point_candidates[i].point_distance_;
121 return k_indices.size();
124template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
130 return (
nearestKSearch(search_point, k, k_indices, k_sqr_distances));
133template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
136 const PointT& p_q,
index_t& result_index,
float& sqr_distance)
140 "Invalid (NaN, Inf) point coordinates given to nearestKSearch!");
143 key.
x = key.
y = key.
z = 0;
146 p_q, this->
root_node_, key, 1, result_index, sqr_distance);
151template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
161template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
167 std::vector<float>& k_sqr_distances,
171 "Invalid (NaN, Inf) point coordinates given to nearestKSearch!");
173 key.
x = key.
y = key.
z = 0;
176 k_sqr_distances.clear();
187 return k_indices.size();
190template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
196 std::vector<float>& k_sqr_distances,
201 return (
radiusSearch(search_point, radius, k_indices, k_sqr_distances, max_nn));
204template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
207 const Eigen::Vector3f& min_pt,
208 const Eigen::Vector3f& max_pt,
213 key.
x = key.
y = key.
z = 0;
219 return k_indices.size();
222template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
231 const double squared_search_radius,
232 std::vector<prioPointQueueEntry>& point_candidates)
const
234 std::vector<prioBranchQueueEntry> search_heap;
235 search_heap.resize(8);
239 double smallest_squared_dist = squared_search_radius;
245 for (
unsigned char child_idx = 0; child_idx < 8; child_idx++) {
249 search_heap[child_idx].key.x = (key.
x << 1) + (!!(child_idx & (1 << 2)));
250 search_heap[child_idx].key.y = (key.
y << 1) + (!!(child_idx & (1 << 1)));
251 search_heap[child_idx].key.z = (key.
z << 1) + (!!(child_idx & (1 << 0)));
255 search_heap[child_idx].key, tree_depth, voxel_center);
259 search_heap[child_idx].point_distance =
pointSquaredDist(voxel_center, point);
262 search_heap[child_idx].point_distance = std::numeric_limits<float>::infinity();
266 std::sort(search_heap.begin(), search_heap.end());
271 while ((!search_heap.empty()) &&
272 (search_heap.back().point_distance <
273 smallest_squared_dist + voxelSquaredDiameter / 4.0 +
274 sqrt(smallest_squared_dist * voxelSquaredDiameter) - this->
epsilon_)) {
278 child_node = search_heap.back().node;
279 new_key = search_heap.back().key;
283 smallest_squared_dist =
289 smallest_squared_dist,
296 const auto* child_leaf =
static_cast<const LeafNode*
>(child_node);
299 (*child_leaf)->getPointIndices(decoded_point_vector);
302 for (
const auto& point_index : decoded_point_vector) {
304 const PointT& candidate_point = this->getPointByIndex(point_index);
307 float squared_dist = pointSquaredDist(candidate_point, point);
309 const auto insert_into_queue = [&] {
310 point_candidates.emplace(
311 std::upper_bound(point_candidates.begin(),
312 point_candidates.end(),
314 [](
float dist,
const prioPointQueueEntry& ent) {
315 return dist < ent.point_distance_;
320 if (point_candidates.size() <
K) {
323 else if (point_candidates.back().point_distance_ > squared_dist) {
324 point_candidates.pop_back();
329 if (point_candidates.size() ==
K)
330 smallest_squared_dist = point_candidates.back().point_distance_;
333 search_heap.pop_back();
336 return (smallest_squared_dist);
339template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
343 const double radiusSquared,
348 std::vector<float>& k_sqr_distances,
352 for (
unsigned char child_idx = 0; child_idx < 8; child_idx++) {
363 new_key.
x = (key.
x << 1) + (!!(child_idx & (1 << 2)));
364 new_key.
y = (key.
y << 1) + (!!(child_idx & (1 << 1)));
365 new_key.
z = (key.
z << 1) + (!!(child_idx & (1 << 0)));
369 Eigen::Vector3f min_pt, max_pt;
372 if (point.x < min_pt.x())
373 squared_dist += std::pow(point.x - min_pt.x(), 2);
374 else if (point.x > max_pt.x())
375 squared_dist += std::pow(point.x - max_pt.x(), 2);
376 if (point.y < min_pt.y())
377 squared_dist += std::pow(point.y - min_pt.y(), 2);
378 else if (point.y > max_pt.y())
379 squared_dist += std::pow(point.y - max_pt.y(), 2);
380 if (point.z < min_pt.z())
381 squared_dist += std::pow(point.z - min_pt.z(), 2);
382 else if (point.z > max_pt.z())
383 squared_dist += std::pow(point.z - max_pt.z(), 2);
384 if (squared_dist < (radiusSquared + this->
epsilon_)) {
395 if (max_nn != 0 && k_indices.size() == max_nn)
400 const auto* child_leaf =
static_cast<const LeafNode*
>(child_node);
403 for (
const auto& index : (*child_leaf)->getPointIndicesVector()) {
404 const PointT& candidate_point = this->getPointByIndex(index);
407 squared_dist = pointSquaredDist(candidate_point, point);
410 if (squared_dist > radiusSquared)
414 k_indices.push_back(index);
415 k_sqr_distances.push_back(squared_dist);
417 if (max_nn != 0 && k_indices.size() == max_nn)
425template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
441 double min_voxel_center_distance = std::numeric_limits<double>::max();
443 unsigned char min_child_idx = 0xFF;
446 for (
unsigned char child_idx = 0; child_idx < 8; child_idx++) {
451 double voxelPointDist;
453 new_key.
x = (key.
x << 1) + (!!(child_idx & (1 << 2)));
454 new_key.
y = (key.
y << 1) + (!!(child_idx & (1 << 1)));
455 new_key.
z = (key.
z << 1) + (!!(child_idx & (1 << 0)));
463 if (voxelPointDist >= min_voxel_center_distance)
466 min_voxel_center_distance = voxelPointDist;
467 min_child_idx = child_idx;
468 minChildKey = new_key;
472 assert(min_child_idx < 8);
474 child_node = this->getBranchChildPtr(*node, min_child_idx);
489 const auto* child_leaf =
static_cast<const LeafNode*
>(child_node);
491 float smallest_squared_dist = std::numeric_limits<float>::max();
494 (**child_leaf).getPointIndices(decoded_point_vector);
497 for (
const auto& index : decoded_point_vector) {
498 const PointT& candidate_point = this->getPointByIndex(index);
504 if (squared_dist >= smallest_squared_dist)
507 result_index = index;
508 smallest_squared_dist = squared_dist;
509 sqr_distance = squared_dist;
514template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
517 const PointT& point_a,
const PointT& point_b)
const
519 return (point_a.getVector3fMap() - point_b.getVector3fMap()).squaredNorm();
522template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
525 const Eigen::Vector3f& min_pt,
526 const Eigen::Vector3f& max_pt,
533 for (
unsigned char child_idx = 0; child_idx < 8; child_idx++) {
543 new_key.
x = (key.
x << 1) + (!!(child_idx & (1 << 2)));
544 new_key.
y = (key.
y << 1) + (!!(child_idx & (1 << 1)));
545 new_key.
z = (key.
z << 1) + (!!(child_idx & (1 << 0)));
548 Eigen::Vector3f lower_voxel_corner;
549 Eigen::Vector3f upper_voxel_corner;
552 new_key, tree_depth, lower_voxel_corner, upper_voxel_corner);
556 if ((lower_voxel_corner(0) <= max_pt(0)) && (min_pt(0) <= upper_voxel_corner(0)) &&
557 (lower_voxel_corner(1) <= max_pt(1)) && (min_pt(1) <= upper_voxel_corner(1)) &&
558 (lower_voxel_corner(2) <= max_pt(2)) && (min_pt(2) <= upper_voxel_corner(2))) {
573 const auto* child_leaf =
static_cast<const LeafNode*
>(child_node);
576 (**child_leaf).getPointIndices(decoded_point_vector);
579 for (
const auto& index : decoded_point_vector) {
584 ((candidate_point.x >= min_pt(0)) && (candidate_point.x <= max_pt(0)) &&
585 (candidate_point.y >= min_pt(1)) && (candidate_point.y <= max_pt(1)) &&
586 (candidate_point.z >= min_pt(2)) && (candidate_point.z <= max_pt(2)));
590 k_indices.push_back(index);
597template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
601 Eigen::Vector3f direction,
606 key.
x = key.
y = key.
z = 0;
608 voxel_center_list.clear();
613 double min_x, min_y, min_z, max_x, max_y, max_z;
617 if (std::max(std::max(min_x, min_y), min_z) < std::min(std::min(max_x, max_y), max_z))
633template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
637 Eigen::Vector3f direction,
642 key.
x = key.
y = key.
z = 0;
648 double min_x, min_y, min_z, max_x, max_y, max_z;
652 if (std::max(std::max(min_x, min_y), min_z) < std::min(std::min(max_x, max_y), max_z))
667template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
682 if (max_x < 0.0 || max_y < 0.0 || max_z < 0.0)
691 voxel_center_list.push_back(newPoint);
700 double mid_x = 0.5 * (min_x + max_x);
701 double mid_y = 0.5 * (min_y + max_y);
702 double mid_z = 0.5 * (min_z + max_z);
708 unsigned char child_idx;
713 child_idx =
static_cast<unsigned char>(curr_node ^ a);
722 child_key.
x = (key.
x << 1) | (!!(child_idx & (1 << 2)));
723 child_key.
y = (key.
y << 1) | (!!(child_idx & (1 << 1)));
724 child_key.
z = (key.
z << 1) | (!!(child_idx & (1 << 0)));
859 }
while ((curr_node < 8) && (max_voxel_count <= 0 || voxel_count < max_voxel_count));
860 return (voxel_count);
863template <
typename Po
intT,
typename LeafContainerT,
typename BranchContainerT>
878 if (max_x < 0.0 || max_y < 0.0 || max_z < 0.0)
883 const auto* leaf =
static_cast<const LeafNode*
>(node);
886 (*leaf)->getPointIndices(k_indices);
895 double mid_x = 0.5 * (min_x + max_x);
896 double mid_y = 0.5 * (min_y + max_y);
897 double mid_z = 0.5 * (min_z + max_z);
903 unsigned char child_idx;
907 child_idx =
static_cast<unsigned char>(curr_node ^ a);
915 child_key.
x = (key.
x << 1) | (!!(child_idx & (1 << 2)));
916 child_key.
y = (key.
y << 1) | (!!(child_idx & (1 << 1)));
917 child_key.
z = (key.
z << 1) | (!!(child_idx & (1 << 0)));
1051 }
while ((curr_node < 8) && (max_voxel_count <= 0 || voxel_count < max_voxel_count));
1053 return (voxel_count);
1059#define PCL_INSTANTIATE_OctreePointCloudSearch(T) \
1060 template class PCL_EXPORTS pcl::octree::OctreePointCloudSearch<T>;
std::size_t leaf_count_
Amount of leaf nodes.
BranchNode * root_node_
Pointer to root branch node of octree.
bool branchHasChild(const BranchNode &branch_arg, unsigned char child_idx_arg) const
Check if branch is pointing to a particular child node.
LeafContainerT * findLeaf(uindex_t idx_x_arg, uindex_t idx_y_arg, uindex_t idx_z_arg) const
Find leaf node at (idx_x_arg, idx_y_arg, idx_z_arg).
OctreeNode * getBranchChildPtr(const BranchNode &branch_arg, unsigned char child_idx_arg) const
Retrieve a child node pointer for child node at child_idx.
Abstract octree node class
virtual node_type_t getNodeType() const =0
Pure virtual method for retrieving the type of octree node (branch or leaf).
const PointT & getPointByIndex(uindex_t index_arg) const
void genLeafNodeCenterFromOctreeKey(const OctreeKey &key_arg, PointT &point_arg) const
void genOctreeKeyforPoint(const PointT &point_arg, OctreeKey &key_arg) const
void genVoxelCenterFromOctreeKey(const OctreeKey &key_arg, uindex_t tree_depth_arg, PointT &point_arg) const
double getVoxelSquaredDiameter(uindex_t tree_depth_arg) const
void genVoxelBoundsFromOctreeKey(const OctreeKey &key_arg, uindex_t tree_depth_arg, Eigen::Vector3f &min_pt, Eigen::Vector3f &max_pt) const
typename OctreeT::LeafNode LeafNode
double getKNearestNeighborRecursive(const PointT &point, uindex_t K, const BranchNode *node, const OctreeKey &key, uindex_t tree_depth, const double squared_search_radius, std::vector< prioPointQueueEntry > &point_candidates) const
uindex_t getIntersectedVoxelIndices(Eigen::Vector3f origin, Eigen::Vector3f direction, Indices &k_indices, uindex_t max_voxel_count=0) const
Get indices of all voxels that are intersected by a ray (origin, direction).
int getNextIntersectedNode(double x, double y, double z, int a, int b, int c) const
Get the next visited node given the current node upper bounding box corner.
uindex_t getIntersectedVoxelIndicesRecursive(double min_x, double min_y, double min_z, double max_x, double max_y, double max_z, unsigned char a, const OctreeNode *node, const OctreeKey &key, Indices &k_indices, uindex_t max_voxel_count) const
Recursively search the tree for all intersected leaf nodes and return a vector of indices.
uindex_t getIntersectedVoxelCenters(Eigen::Vector3f origin, Eigen::Vector3f direction, AlignedPointTVector &voxel_center_list, uindex_t max_voxel_count=0) const
Get a PointT vector of centers of all voxels that intersected by a ray (origin, direction).
bool voxelSearch(const PointT &point, Indices &point_idx_data)
Search for neighbors within a voxel at given point.
std::vector< PointT, Eigen::aligned_allocator< PointT > > AlignedPointTVector
uindex_t boxSearch(const Eigen::Vector3f &min_pt, const Eigen::Vector3f &max_pt, Indices &k_indices) const
Search for points within rectangular search area Points exactly on the edges of the search rectangle ...
uindex_t getIntersectedVoxelCentersRecursive(double min_x, double min_y, double min_z, double max_x, double max_y, double max_z, unsigned char a, const OctreeNode *node, const OctreeKey &key, AlignedPointTVector &voxel_center_list, uindex_t max_voxel_count) const
Recursively search the tree for all intersected leaf nodes and return a vector of voxel centers.
void approxNearestSearchRecursive(const PointT &point, const BranchNode *node, const OctreeKey &key, uindex_t tree_depth, index_t &result_index, float &sqr_distance)
Recursive search method that explores the octree and finds the approximate nearest neighbor.
uindex_t nearestKSearch(const PointCloud &cloud, uindex_t index, uindex_t k, Indices &k_indices, std::vector< float > &k_sqr_distances)
Search for k-nearest neighbors at the query point.
void boxSearchRecursive(const Eigen::Vector3f &min_pt, const Eigen::Vector3f &max_pt, const BranchNode *node, const OctreeKey &key, uindex_t tree_depth, Indices &k_indices) const
Recursive search method that explores the octree and finds points within a rectangular search area.
void initIntersectedVoxel(Eigen::Vector3f &origin, Eigen::Vector3f &direction, double &min_x, double &min_y, double &min_z, double &max_x, double &max_y, double &max_z, unsigned char &a) const
Initialize raytracing algorithm.
float pointSquaredDist(const PointT &point_a, const PointT &point_b) const
void getNeighborsWithinRadiusRecursive(const PointT &point, const double radiusSquared, const BranchNode *node, const OctreeKey &key, uindex_t tree_depth, Indices &k_indices, std::vector< float > &k_sqr_distances, uindex_t max_nn) const
Recursive search method that explores the octree and finds neighbors within a given radius.
typename OctreeT::BranchNode BranchNode
void approxNearestSearch(const PointCloud &cloud, uindex_t query_index, index_t &result_index, float &sqr_distance)
Search for approx.
int getFirstIntersectedNode(double min_x, double min_y, double min_z, double mid_x, double mid_y, double mid_z) const
Find first child node ray will enter.
uindex_t radiusSearch(const PointCloud &cloud, uindex_t index, double radius, Indices &k_indices, std::vector< float > &k_sqr_distances, index_t max_nn=0)
Search for all neighbors of query point that are within a given radius.
bool isFinite(const PointT &pt)
Tests if the 3D components of a point are all finite param[in] pt point to be tested return true if f...
detail::int_type_t< detail::index_type_size, detail::index_type_signed > index_t
Type used for an index in PCL.
IndicesAllocator<> Indices
Type used for indices in PCL.
detail::int_type_t< detail::index_type_size, false > uindex_t
Type used for an unsigned index in PCL.