From 8d2a13ecbbf5012d7bd243457a23e79c4081b8a2 Mon Sep 17 00:00:00 2001 From: Jaewook Lee Date: Wed, 15 Nov 2023 16:57:37 +0100 Subject: [PATCH 1/3] look for valid nanoflann in system --- CMakeLists.txt | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 15a4a36..ebba5c3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.0) +cmake_minimum_required(VERSION 3.12) project(napf VERSION 0.0.0 LANGUAGES CXX) if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) @@ -22,18 +22,27 @@ set(TARGETS_EXPORT_NAME "${PROJECT_NAME}Targets") set(namespace "${PROJECT_NAME}::") # sources -set(CXX_HEADERS src/napf.hpp third_party/nanoflann.hpp) +set(CXX_HEADERS src/napf.hpp) # Interface Library, since it's header only lib add_library(napf INTERFACE) add_library(napf::napf ALIAS napf) -target_include_directories(napf - INTERFACE - $ - $ - $) + +# try to use installed nanoflann. +# else, include the one in the third_party +find_package(nanoflann QUIET) +if(nanoflann_FOUND AND nanoflann_VERSION VERSION_GREATER_EQUAL "1.5.0") + message("nanoflann found - napf will link to nanoflann found in system.") + target_link_libraries(napf INTERFACE nanoflann::nanoflann) +else() + set(CXX_HEADERS ${CXX_HEADERS} third_party/nanoflann.hpp) + target_include_directories(napf + INTERFACE + $ + $ + $) +endif() target_compile_features(napf INTERFACE cxx_std_11) -#target_compile_options(napf PRIVATE -O3) # everyone can enjoy optimization if(BUILD_FORTRAN_MODULE) message("*** building additional fortran module ***") From 7fcba64613affeafc8761665cecd19a98d57dbda Mon Sep 17 00:00:00 2001 From: Jaewook Lee Date: Tue, 2 Jan 2024 21:34:33 +0100 Subject: [PATCH 2/3] set basic include path --- CMakeLists.txt | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ebba5c3..4f81784 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -28,6 +28,12 @@ set(CXX_HEADERS src/napf.hpp) add_library(napf INTERFACE) add_library(napf::napf ALIAS napf) +# basic path +target_include_directories(napf + INTERFACE + $ + $) + # try to use installed nanoflann. # else, include the one in the third_party find_package(nanoflann QUIET) @@ -38,9 +44,7 @@ else() set(CXX_HEADERS ${CXX_HEADERS} third_party/nanoflann.hpp) target_include_directories(napf INTERFACE - $ - $ - $) + $) endif() target_compile_features(napf INTERFACE cxx_std_11) From 941705cd4d3b4069116209bba8b90debe77a2ec0 Mon Sep 17 00:00:00 2001 From: Jaewook Lee Date: Tue, 2 Jan 2024 21:36:17 +0100 Subject: [PATCH 3/3] update nanoflann --- third_party/nanoflann | 2 +- third_party/nanoflann.hpp | 158 ++++++++++++++++++++++++++++++++------ 2 files changed, 136 insertions(+), 24 deletions(-) diff --git a/third_party/nanoflann b/third_party/nanoflann index f1e47f0..0b9fc58 160000 --- a/third_party/nanoflann +++ b/third_party/nanoflann @@ -1 +1 @@ -Subproject commit f1e47f0dbae72c11b008fd0d5b1191d9ba1f29bb +Subproject commit 0b9fc585e3425e821c834ff8e4e584f6f8007f63 diff --git a/third_party/nanoflann.hpp b/third_party/nanoflann.hpp index c2c80fd..b6d3315 100644 --- a/third_party/nanoflann.hpp +++ b/third_party/nanoflann.hpp @@ -3,7 +3,7 @@ * * Copyright 2008-2009 Marius Muja (mariusm@cs.ubc.ca). All rights reserved. * Copyright 2008-2009 David G. Lowe (lowe@cs.ubc.ca). All rights reserved. - * Copyright 2011-2022 Jose Luis Blanco (joseluisblancoc@gmail.com). + * Copyright 2011-2023 Jose Luis Blanco (joseluisblancoc@gmail.com). * All rights reserved. * * THE BSD LICENSE @@ -50,6 +50,7 @@ #include #include // for abs() #include // for abs() +#include #include // std::reference_wrapper #include #include @@ -60,7 +61,7 @@ #include /** Library version: 0xMmP (M=Major,m=minor,P=patch) */ -#define NANOFLANN_VERSION 0x150 +#define NANOFLANN_VERSION 0x151 // Avoid conflicting declaration of min/max macros in Windows headers #if !defined(NOMINMAX) && \ @@ -158,6 +159,8 @@ inline typename std::enable_if::value, void>::type /** @addtogroup result_sets_grp Result set classes * @{ */ + +/** Result set for KNN searches (N-closest neighbors) */ template < typename _DistanceType, typename _IndexType = size_t, typename _CountType = size_t> @@ -190,8 +193,93 @@ class KNNResultSet } CountType size() const { return count; } + bool empty() const { return count == 0; } + bool full() const { return count == capacity; } + + /** + * Called during search to add an element matching the criteria. + * @return true if the search should be continued, false if the results are + * sufficient + */ + bool addPoint(DistanceType dist, IndexType index) + { + CountType i; + for (i = count; i > 0; --i) + { + /** If defined and two points have the same distance, the one with + * the lowest-index will be returned first. */ +#ifdef NANOFLANN_FIRST_MATCH + if ((dists[i - 1] > dist) || + ((dist == dists[i - 1]) && (indices[i - 1] > index))) + { +#else + if (dists[i - 1] > dist) + { +#endif + if (i < capacity) + { + dists[i] = dists[i - 1]; + indices[i] = indices[i - 1]; + } + } + else + break; + } + if (i < capacity) + { + dists[i] = dist; + indices[i] = index; + } + if (count < capacity) count++; + + // tell caller that the search shall continue + return true; + } + + DistanceType worstDist() const { return dists[capacity - 1]; } +}; + +/** Result set for RKNN searches (N-closest neighbors with a maximum radius) */ +template < + typename _DistanceType, typename _IndexType = size_t, + typename _CountType = size_t> +class RKNNResultSet +{ + public: + using DistanceType = _DistanceType; + using IndexType = _IndexType; + using CountType = _CountType; - bool full() const { return count == capacity; } + private: + IndexType* indices; + DistanceType* dists; + CountType capacity; + CountType count; + DistanceType maximumSearchDistanceSquared; + + public: + explicit RKNNResultSet( + CountType capacity_, DistanceType maximumSearchDistanceSquared_) + : indices(nullptr), + dists(nullptr), + capacity(capacity_), + count(0), + maximumSearchDistanceSquared(maximumSearchDistanceSquared_) + { + } + + void init(IndexType* indices_, DistanceType* dists_) + { + indices = indices_; + dists = dists_; + count = 0; + if (capacity) + dists[capacity - 1] = maximumSearchDistanceSquared; + } + + CountType size() const { return count; } + bool empty() const { return count == 0; } + bool full() const { return count == capacity; } /** * Called during search to add an element matching the criteria. @@ -744,7 +832,7 @@ struct SearchParameters */ class PooledAllocator { - static constexpr size_t WORDSIZE = 16; + static constexpr size_t WORDSIZE = 16; //WORDSIZE must >= 8 static constexpr size_t BLOCKSIZE = 8192; /* We maintain memory alignment to word boundaries by requiring that all @@ -753,9 +841,7 @@ class PooledAllocator /* Minimum number of bytes requested at a time from the system. Must be * multiple of WORDSIZE. */ - using Offset = uint32_t; - using Size = uint32_t; - using Dimension = int32_t; + using Size = size_t; Size remaining_ = 0; //!< Number of bytes left in current block of storage void* base_ = nullptr; //!< Pointer to base of current block of storage @@ -817,9 +903,9 @@ class PooledAllocator /* Allocate new storage. */ const Size blocksize = - (size + sizeof(void*) + (WORDSIZE - 1) > BLOCKSIZE) - ? size + sizeof(void*) + (WORDSIZE - 1) - : BLOCKSIZE; + size > BLOCKSIZE + ? size + WORDSIZE + : BLOCKSIZE + WORDSIZE; // use the standard C malloc to allocate memory void* m = ::malloc(blocksize); @@ -833,12 +919,8 @@ class PooledAllocator static_cast(m)[0] = base_; base_ = m; - Size shift = 0; - // int size_t = (WORDSIZE - ( (((size_t)m) + sizeof(void*)) & - // (WORDSIZE-1))) & (WORDSIZE-1); - - remaining_ = blocksize - sizeof(void*) - shift; - loc_ = (static_cast(m) + sizeof(void*) + shift); + remaining_ = blocksize - WORDSIZE; + loc_ = static_cast(m) + WORDSIZE; } void* rloc = loc_; loc_ = static_cast(loc_) + size; @@ -1220,25 +1302,26 @@ class KDTreeBaseClass } ElementType max_spread = -1; cutfeat = 0; + ElementType min_elem = 0, max_elem = 0; for (Dimension i = 0; i < dims; ++i) { ElementType span = bbox[i].high - bbox[i].low; if (span > (1 - EPS) * max_span) { - ElementType min_elem, max_elem; - computeMinMax(obj, ind, count, i, min_elem, max_elem); - ElementType spread = max_elem - min_elem; + ElementType min_elem_, max_elem_; + computeMinMax(obj, ind, count, i, min_elem_, max_elem_); + ElementType spread = max_elem_ - min_elem_; if (spread > max_spread) { cutfeat = i; max_spread = spread; + min_elem = min_elem_; + max_elem = max_elem_; } } } // split in the middle DistanceType split_val = (bbox[cutfeat].low + bbox[cutfeat].high) / 2; - ElementType min_elem, max_elem; - computeMinMax(obj, ind, count, cutfeat, min_elem, max_elem); if (split_val < min_elem) cutval = min_elem; @@ -1680,6 +1763,34 @@ class KDTreeSingleIndexAdaptor return resultSet.size(); } + /** + * Find the first N neighbors to \a query_point[0:dim-1] within a maximum + * radius. The output is given as a vector of pairs, of which the first + * element is a point index and the second the corresponding distance. + * Previous contents of \a IndicesDists are cleared. + * + * \sa radiusSearch, findNeighbors + * \return Number `N` of valid points in the result set. + * + * \note If L2 norms are used, all returned distances are actually squared + * distances. + * + * \note Only the first `N` entries in `out_indices` and `out_distances` + * will be valid. Return is less than `num_closest` only if the + * number of elements in the tree is less than `num_closest`. + */ + Size rknnSearch( + const ElementType* query_point, const Size num_closest, + IndexType* out_indices, DistanceType* out_distances, + const DistanceType& radius) const + { + nanoflann::RKNNResultSet resultSet( + num_closest, radius); + resultSet.init(out_indices, out_distances); + findNeighbors(resultSet, query_point); + return resultSet.size(); + } + /** @} */ public: @@ -2060,11 +2171,12 @@ class KDTreeSingleIndexDynamicAdaptor_ */ Size knnSearch( const ElementType* query_point, const Size num_closest, - IndexType* out_indices, DistanceType* out_distances) const + IndexType* out_indices, DistanceType* out_distances, + const SearchParameters& searchParams = {}) const { nanoflann::KNNResultSet resultSet(num_closest); resultSet.init(out_indices, out_distances); - findNeighbors(resultSet, query_point); + findNeighbors(resultSet, query_point, searchParams); return resultSet.size(); }