Repository: drprojects/point_geometric_features (language: python)

Issue:
pgeof as a pip package
Packaging `pgeof` as a pip package, free from the conda dependency of `eigen3` would be great.

@rjanvier I think you were thinking of working on this at some point. Is this still something you are considering ?
Multi scale features computation
Hi Damien,
I have a little window to try to implement multi scale feature computation we talked about (see #4).
I would like to know how do you see the thing? would you want to rely on an aggregation scheme (I don't think a mean  makes sense for all features) or output the whole multi scale feature map?
Thank in advance,
Romain


Discussion / hints:

Hi @rjanvier, thanks for reaching out :blush: 

The way I see it, there is no need for an aggregation scheme. The trick would be simply to compute the same handcrafted features we already do. Except that we could use different subsets of the provided neighbors. 

Assuming the neighbors are passed in order of increasing distance, we can easily compute the features for a series of increasing scales. These could either be defined as K-NN scales (eg $K \in (10, 25, 50, ... ) $) or radius scales (eg $R \in ( 0.5, 0.1, 0.2, ... ) $). 

Do you see what I mean ? Were you thinking of anything different ? 
Hi Damien, thank you for your quick answer.

Yes, I think I follow you and I already have this "trick" in mind since the input requirements of pgeof easily allow this behavior. 

My question is more about the output format / downstream usage of the multiscale features. If I follow your last post, we will output a feature map for each K / R selected by user so there is no aggregation. Usage of multiple scales or a putative aggregation is left to the responsability of the "caller" of pgeof am I right?
Absolutely ! I am thinking that it is up to the user to decide what to do with the resulting features. This should not bring too much change to the code base, only need to take a list of K / R as input.

Thanks a lot for looking into this ! 

Fix patch:
diff --git a/.clang-format b/.clang-format
new file mode 100644
index 0000000..aa6a093
--- /dev/null
+++ b/.clang-format
@@ -0,0 +1,83 @@
+Language:        Cpp
+BasedOnStyle: Google
+# ---
+#AccessModifierOffset: -4
+AlignAfterOpenBracket: AlwaysBreak # Values: Align, DontAlign, AlwaysBreak
+AlignConsecutiveAssignments: true
+AlignConsecutiveDeclarations: true
+#AlignEscapedNewlinesLeft: true
+#AlignOperands:   false
+AlignTrailingComments: false # Should be off, causes many dummy problems!!
+#AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortBlocksOnASingleLine: true
+#AllowShortCaseLabelsOnASingleLine: false
+#AllowShortFunctionsOnASingleLine: Empty
+#AllowShortIfStatementsOnASingleLine: false
+#AllowShortLoopsOnASingleLine: false
+#AlwaysBreakAfterDefinitionReturnType: None
+#AlwaysBreakAfterReturnType: None
+#AlwaysBreakBeforeMultilineStrings: true
+#AlwaysBreakTemplateDeclarations: true
+#BinPackArguments: false
+#BinPackParameters: false
+#BraceWrapping:
+  #AfterClass:      false
+  #AfterControlStatement: false
+  #AfterEnum:       false
+  #AfterFunction:   false
+  #AfterNamespace:  false
+  #AfterObjCDeclaration: false
+  #AfterStruct:     false
+  #AfterUnion:      false
+  #BeforeCatch:     false
+  #BeforeElse:      true
+  #IndentBraces:    false
+#BreakBeforeBinaryOperators: None
+BreakBeforeBraces: Allman
+#BreakBeforeTernaryOperators: true
+#BreakConstructorInitializersBeforeComma: false
+ColumnLimit: 120
+#CommentPragmas:  ''
+#ConstructorInitializerAllOnOneLineOrOnePerLine: true
+#ConstructorInitializerIndentWidth: 4
+#ContinuationIndentWidth: 4
+#Cpp11BracedListStyle: true
+#DerivePointerAlignment: false
+#DisableFormat:   false
+#ExperimentalAutoDetectBinPacking: false
+##FixNamespaceComments: true # Not applicable in 3.8
+#ForEachMacros:   [ foreach, Q_FOREACH, BOOST_FOREACH ]
+#IncludeCategories:
+  #- Regex:           '.*'
+    #Priority:        1
+IndentCaseLabels: true
+IndentWidth:     4
+IndentWrappedFunctionNames: true
+#KeepEmptyLinesAtTheStartOfBlocks: true
+#MacroBlockBegin: ''
+#MacroBlockEnd:   ''
+MaxEmptyLinesToKeep: 1
+NamespaceIndentation: None
+#PenaltyBreakBeforeFirstCallParameter: 19
+#PenaltyBreakComment: 300
+#PenaltyBreakFirstLessLess: 120
+#PenaltyBreakString: 1000
+#PenaltyExcessCharacter: 1000000
+#PenaltyReturnTypeOnItsOwnLine: 200
+DerivePointerAlignment: false
+#PointerAlignment: Left
+ReflowComments:  true # Should be true, otherwise clang-format doesn't touch comments
+SortIncludes:    true
+#SpaceAfterCStyleCast: false
+SpaceBeforeAssignmentOperators: true
+#SpaceBeforeParens: ControlStatements
+#SpaceInEmptyParentheses: false
+#SpacesBeforeTrailingComments: 2
+#SpacesInAngles:  false
+#SpacesInContainerLiterals: true
+#SpacesInCStyleCastParentheses: false
+#SpacesInParentheses: false
+#SpacesInSquareBrackets: false
+Standard:        Cpp11
+TabWidth:        4
+UseTab:          Never # Available options are Never, Always, ForIndentation
diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml
index e972cc2..6cf119a 100644
--- a/.github/workflows/wheels.yml
+++ b/.github/workflows/wheels.yml
@@ -15,21 +15,16 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        os: [ubuntu-20.04, windows-2019]
+        os: [ubuntu-22.04, windows-2022, macos-12]
+      fail-fast: false
 
-    env:
-      CIBW_ENVIRONMENT:  ${{ matrix.os == 'windows-2019' && 'EIGEN_LIB_PATH=.' || 'EIGEN_LIB_PATH=/project' }}
     steps:
-      - uses: actions/checkout@v4
-
       - uses: actions/checkout@v4
         with:
-          repository: eigen-mirror/eigen
-          ref: "3.4.0"
-          path: eigen3
-
+          submodules: True
+        
       # Used to host cibuildwheel
-      - uses: actions/setup-python@v3
+      - uses: actions/setup-python@v4
 
       - name: Install cibuildwheel
         run: python -m pip install cibuildwheel==2.15
@@ -49,9 +44,11 @@ jobs:
     # Only publish to PyPI when a commit is tagged
     if: github.event_name == 'push' && startsWith(github.event.ref, 'refs/tags')
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
+        with:
+          submodules: True
 
-      - uses: actions/setup-python@v3
+      - uses: actions/setup-python@v4
 
       - uses: actions/download-artifact@v3
         with:
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..26d6e71
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,10 @@
+[submodule "eigen"]
+	path = third_party/eigen
+	url = https://gitlab.com/libeigen/eigen
+	branch = 3.4
+[submodule "nanoflann"]
+	path = third_party/nanoflann
+	url = https://github.com/jlblancoc/nanoflann
+[submodule "taskflow"]
+	path = third_party/taskflow
+	url = https://github.com/taskflow/taskflow
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..b8ee37a
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,33 @@
+cmake_minimum_required(VERSION 3.15...3.27)
+project(pgeof)
+find_package(Python 3.8 COMPONENTS Interpreter Development.Module REQUIRED)
+
+if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
+  set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build." FORCE)
+  set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo")
+endif()
+
+# Detect the installed nanobind package and import it into CMake
+execute_process(
+  COMMAND "${Python_EXECUTABLE}" -m nanobind --cmake_dir
+  OUTPUT_STRIP_TRAILING_WHITESPACE OUTPUT_VARIABLE NB_DIR)
+list(APPEND CMAKE_PREFIX_PATH "${NB_DIR}")
+find_package(nanobind CONFIG REQUIRED)
+
+nanobind_add_module(pgeof_ext NOMINSIZE STABLE_ABI LTO src/pgeof_ext.cpp)
+
+# for nanobind 2
+# nanobind_add_stub(
+#   pgeof_ext_stub
+#   MODULE pgeof_ext
+#   OUTPUT pgeof_ext.pyi
+#   PYTHON_PATH $<TARGET_FILE_DIR:pgeof_ext>
+#   DEPENDS pgeof_ext
+# )
+
+# All lib are headeer only. 
+# it's faster to include like this than using exported targets 
+# (i.e add_subdirectories(...))
+target_include_directories(pgeof_ext PRIVATE "include" "third_party/eigen" "third_party/nanoflann/include" "third_party/taskflow") 
+
+install(TARGETS pgeof_ext LIBRARY DESTINATION pgeof)
\ No newline at end of file
diff --git a/MANIFEST.in b/MANIFEST.in
deleted file mode 100644
index 19e15b6..0000000
--- a/MANIFEST.in
+++ /dev/null
@@ -1,1 +0,0 @@
-graft include
\ No newline at end of file
diff --git a/README.md b/README.md
index a53a93d..f0d85a4 100644
--- a/README.md
+++ b/README.md
@@ -2,10 +2,9 @@
 
 # Point Geometric Features
 
-![Linux](https://img.shields.io/badge/Linux-FCC624?style=for-the-badge&logo=linux&logoColor=black)
-[![python](https://img.shields.io/badge/-Python_3.8-blue?logo=python&logoColor=white)](https://github.com/pre-commit/pre-commit)
+[![python](https://img.shields.io/badge/-Python_3.9_%7C_3.10_%7C_3.11_%7C_3.12-blue?logo=python&logoColor=white)](#)
 ![C++](https://img.shields.io/badge/c++-%2300599C.svg?style=for-the-badge&logo=c%2B%2B&logoColor=white)
-[![license](https://img.shields.io/badge/License-MIT-green.svg?labelColor=gray)](https://github.com/ashleve/lightning-hydra-template#license)
+[![license](https://img.shields.io/badge/License-MIT-green.svg?labelColor=gray)](#)
 
 
 </div>
@@ -21,7 +20,7 @@ Python wrapper around C++ helper to compute, for each point in a 3D point cloud,
 - linearity
 - planarity
 - scattering
-- verticality
+- verticality (two formulations)
 - normal_x
 - normal_y
 - normal_z
@@ -32,109 +31,154 @@ Python wrapper around C++ helper to compute, for each point in a 3D point cloud,
 - optimal neighborhood size
 </details>
 
+The wrapper allows to compute feature in multiple fashions (on the fly subset of features _a la_ jakteristics, an array of features or multiscale features...). Moreover, it offers basic interfaces to compute fast K-NN or Radius search on point clouds. 
+The overall code is not intended to be DRY nor generic, it aims at providing efficient as possible implementations for some limited scopes and usages.
 
 ## 🧱 Installation
 
-Pgeof will __soon__ be available as pre compiled package on PyPI for both Linux and Windows OSes.
-
 ```bash
 python -m pip install pgeof 
 ```
 
+or 
+
+```bash
+python -m pip install git+https://github.com/drprojects/point_geometric_features
+```
+
 ### building from sources
 
-Pgeof depends on [Eigen library](https://eigen.tuxfamily.org/) and numpy headers at build time.
-The good version of numpy will be fetched from PyPI automatically by the build system but your are responsible for providing
-the path to the Eigen library you want to use (for example py using `CXXFLAGS` variable on Linux or setting `EIGEN_LIB_PATH`)
+pgeof depends on [Eigen library](https://eigen.tuxfamily.org/), [Taskflow](https://github.com/taskflow/taskflow), 
+[nanoflann](https://github.com/jlblancoc/nanoflann) and [nanobind](https://github.com/wjakob/nanobind).
+
+
+pgeof adhere to [PEP 517](https://peps.python.org/pep-0517/) and use [scikit-build-core](https://github.com/scikit-build/scikit-build-core) as build backend. Build dependencies (nanobind, scikit-build-core...) are fetched at build time. C++ third party libraries are embedded as submodules.
+
 
 ```bash
 # clone project
-git clone https://github.com/drprojects/point_geometric_features.git
+git clone --recurse-submodules https://github.com/drprojects/point_geometric_features.git
 cd point_geometric_features
-
-# set the EIGEN_LIB_PATH if needed
-export EIGEN_LIB_PATH="path_to_eigen_root_dir"
 # build and install the package
 python -m pip install .
 ```
 
-### conda 
+## 🚀 Using Point Geometric Features
 
-The following will install the project in a new `pgeof` conda environment.
+Here we summarize the very basics of `pgeof` usage. 
+Users are invited to use `help(pgeof)` for further details on parameters.
 
-```bash
-# clone project
-git clone https://github.com/drprojects/point_geometric_features.git
-cd point_geometric_features
+At its core `pgeof` provides three functions to compute a set of features given a 3D point cloud and
+some precomputed neighborhoods.
 
-# Installation in a new dedicated `pgeof` conda environment
-bash install.sh
+```python
+import pgeof
+
+# Compute a set of 11 predefined features per points.
+pgeof.compute_features(
+    xyz, # The point cloud. A numpy array of shape (n, 3).
+    nn, # CSR data structure see below.
+    nn_ptr, # CSR data structure see below.
+    k_min = 1 # Minimum number of neighbors to consider for features computation. 
+    verbose = false # Basic verbose output, for debug purposes.
+)
 ```
 
-You can easily adapt `install.sh` to install the project in an already-existing 
-environment.
-
-## 🚀 Using Point Geometric Features
-
-The `pgeof` function should be used as follows:
+```python
+# Sequence of n scales feature computation.
+pgeof.compute_features_multiscale(
+    ...
+    k_scale # array of neighborhood size
+)
+```
 
 ```python
-from pgeof import pgeof
-
-pgeof(
-    xyz,              # [n_points, 3] float32 2D array - 3D point coordinates
-    nn,               # [num_neighborhoods] uint32 1D array - Flattened neighbor indices. Make sure those are all positive, '-1' indices will either crash or silently compute incorrect features
-    nn_ptr,           # [n_points+1] uint32 1D array - Pointers wrt `nn`. More specifically, the neighbors of point `i` are `nn[nn_ptr[i]:nn_ptr[i + 1]]`
-    k_min=1,          # (optional, default=1) int - Minimum number of neighbors to consider for features computation. If a point has less, it will be given 0 features
-    k_step=-1,        # (optional, default=-1) int - Step size to take when searching for the optimal neighborhood size for each point, following: http://lareg.ensg.eu/labos/matis/pdf/articles_revues/2015/isprs_wjhm_15.pdf. If k_step < 1, pgeof will not search for the optimal neighborhood and features will be computed based on the all available neighbors for each point 
-    k_min_search=10,  # (optional, default=10) int - Minimum neighborhood size at which to start when searching for the optimal neighborhood size for each point. It is advised to use a value of 10 or higher, for geometric features robustness
-    verbose=False)    # (optional, default=False) bool - Whether computation progress should be printed out
-
-# Print details on how pgeof works and expected input parameters
-print(help(pgeof))
+# Feature computation with optimal neighborhood selection as exposed in Weinmann et al., 2015. 
+# return a set of 12 features per points (11 + the optimal neighborhood size)
+pgeof.compute_features_optimal(
+    ...
+    k_min = 1, # Minimum number of neighbors to consider for features computation.
+    k_step = 1, # Step size to take when searching for the optimal neighborhood
+    k_min_search = 1, # Minimum neighborhood size at which to start when searching for the optimal neighborhood size for each point. it should be >= to k_min parameter.
+)
 ```
 
-👇 You may check out the provided `demo.py` script to get started.
+⚠️ Please note that for theses three functions the **neighbors are expected in CSR format**. 
+This allows expressing neighborhoods of varying sizes with dense arrays (eg the output of a 
+radius search).
 
-```bash
-python demo.py
-```
+We provide very tiny and specialized k-NN / radius-NN search routines. 
+They rely on `nanoflann` C++ library and they should be faster and lighter than `scipy` and `sklearn` alternatives.
 
-⚠️ Please note the **neighbors are expected in CSR format**. This allows 
-expressing neighborhoods of varying sizes with dense arrays (eg the output of a 
-radius search). Here are examples of how to easily convert typical k-NN or 
-radius-NN neighborhoods to CSR format.
+Here are some examples of how to easily compute and convert typical k-NN or radius-NN neighborhoods to CSR format (`nn` and `nn_ptr` are two flat `uint32` arrays):
 
 ```python
-from sklearn.neighbors import NearestNeighbors
+import pgeof
 import numpy as np
 
 # Generate a random synthetic point cloud and k-nearest neighbors
 num_points = 10000
 k = 20
-xyz = np.random.rand(num_points, 3)
-kneigh = NearestNeighbors(n_neighbors=k).fit(xyz).kneighbors(xyz)
+xyz = np.random.rand(num_points, 3).astype("float32")
+knn, _ = pgeof.knn_search(xyz, xyz, k)
 
 # Converting k-nearest neighbors to CSR format
 nn_ptr = np.arange(num_points + 1) * k
-nn = kneigh[1].flatten()
+nn = knn.flatten()
+
+# You may need to convert nn/nn_ptr to uint32 arrays
+nn_ptr = nn_ptr.astype("uint32")
+nn = nn.astype("uint32")
+
+features = pgeof.compute_features(xyz, nn, nn_ptr)
 ```
 
 ```python
-from sklearn.neighbors import NearestNeighbors
+import pgeof
 import numpy as np
 
-# Generate a random synthetic point cloud and radius neighbors
+# Generate a random synthetic point cloud and k-nearest neighbors
 num_points = 10000
-radius = 0.1
-xyz = np.random.rand(num_points, 3)
-rneigh = NearestNeighbors(radius=radius).fit(xyz).radius_neighbors(xyz)
+radius = 0.2
+k = 20
+xyz = np.random.rand(num_points, 3).astype("float32")
+knn, _ = pgeof.radius_search(xyz, xyz, radius, k)
 
 # Converting radius neighbors to CSR format
-nn_ptr = np.r_[0, np.array([x.shape[0] for x in rneigh[1]]).cumsum()]
-nn = np.concatenate(rneigh[1])
+nn_ptr = np.r_[0, (knn >= 0).sum(axis=1).cumsum()]
+nn = knn[knn >= 0]
+# You may need to convert nn/nn_ptr to uint32 arrays
+...
 ```
 
+At last and as a by product we also provide a function to compute a subset of features on the fly. 
+it is inspired by the `Jakteristics` python package (while being less complete but faster).
+the list of feature to compute is given as an array of `EFeatureID`
+
+```python
+import pgeof
+from pgeof import EFeatureID
+import numpy as np
+
+# Generate a random synthetic point cloud and k-nearest neighbors
+num_points = 10000
+radius = 0.2
+k = 20
+xyz = np.random.rand(num_points, 3)
+features = pgeof.compute_features_selected(xyz, radius, k, [EFeatureID.Verticality])
+```
+
+## Known limitations
+
+Some functions only accept `float` scalar types and `uint32` index types and we avoid implicit cast / conversions.
+This could be a limitation in some situations (`double` coordinates point clouds or need for a big indices).
+Some C++ function could be templated / to accept other type without conversion. For now, this feature is not enabled everywhere to reduce compilation time or enhance code readability but please let us know if you need this feature !
+
+Normal are forced to be oriented in the Z half space. 
+
+## Testing
+
+Some basic tests and benchmarks are provided in the tests directory. Test could be run in a clean and reproducible environments via `tox` (`tox run` and `tox run -e bench`).
 
 ## 💳 Credits
 This implementation was largely inspired from [Superpoint Graph](https://github.com/loicland/superpoint_graph). The main modifications here allow: 
@@ -143,6 +187,17 @@ This implementation was largely inspired from [Superpoint Graph](https://github.
 - optimal neighborhood search from this [paper](http://lareg.ensg.eu/labos/matis/pdf/articles_revues/2015/isprs_wjhm_15.pdf)
 - some corrections on geometric features computation
 
+Some heavy refactoring (port to nanobind, test, benchmarks), packaging, speed optimization, feature addition (NN search, on the fly feature computation...) were funded by:
+
+Centre of Wildfire Research of Swansea University (UK) in collaboration with the Research Institute of Biodiversity (CSIC, Spain) and the Department of Mining Exploitation of the University of Oviedo (Spain).
+
+Funding provided by the UK NERC project (NE/T001194/1):
+
+'Advancing 3D Fuel Mapping for Wildfire Behaviour and Risk Mitigation Modelling'
+
+and by the Spanish Knowledge Generation project (PID2021-126790NB-I00):
+
+‘Advancing carbon emission estimations from wildfires applying artificial intelligence to 3D terrestrial point clouds’.
 
 ## License
 
diff --git a/demo.py b/demo.py
deleted file mode 100644
index 35d714b..0000000
--- a/demo.py
+++ /dev/null
@@ -1,54 +0,0 @@
-from pgeof import pgeof
-import numpy as np
-
-# Generate a random synthetic point cloud
-num_points = 10000
-xyz = np.random.rand(num_points, 3)
-
-# Manually generating random neighbors in CSR format
-nn_ptr = np.r_[0, np.random.randint(low=0, high=30, size=num_points).cumsum()]
-nn = np.random.randint(low=0, high=num_points, size=nn_ptr[-1])
-
-# Converting k-nearest neighbors to CSR format
-from sklearn.neighbors import NearestNeighbors
-k = 20
-kneigh = NearestNeighbors(n_neighbors=k).fit(xyz).kneighbors(xyz)
-nn_ptr = np.arange(num_points + 1) * k
-nn = kneigh[1].flatten()
-
-# Converting radius neighbors to CSR format
-from sklearn.neighbors import NearestNeighbors
-radius = 0.1
-rneigh = NearestNeighbors(radius=radius).fit(xyz).radius_neighbors(xyz)
-nn_ptr = np.r_[0, np.array([x.shape[0] for x in rneigh[1]]).cumsum()]
-nn = np.concatenate(rneigh[1])
-
-# Make sure xyz are float32 and nn and nn_ptr are uint32
-xyz = xyz.astype('float32')
-nn_ptr = nn_ptr.astype('uint32')
-nn = nn.astype('uint32')
-
-# Make sure arrays are contiguous (C-order) and not Fortran-order
-xyz = np.ascontiguousarray(xyz)
-nn_ptr = np.ascontiguousarray(nn_ptr)
-nn = np.ascontiguousarray(nn)
-
-# Print details on how pgeof works and expected input parameters
-print(help(pgeof))
-
-# Features have shape [N, 12]:
-#   0 - linearity
-#   1 - planarity
-#   2 - scattering
-#   3 - verticality
-#   4 - normal_x (oriented towards positive z-coordinates)
-#   5 - normal_y (oriented towards positive z-coordinates)
-#   6 - normal_z (oriented towards positive z-coordinates)
-#   7 - length
-#   8 - surface
-#   9 - volume
-#  10 - curvature
-#  11 - optimal neighborhood size (if 'k_step' argument is used)
-geof = pgeof(
-    xyz, nn, nn_ptr, k_min=10, k_step=1, k_min_search=15,
-    verbose=True)
\ No newline at end of file
diff --git a/include/nn_search.hpp b/include/nn_search.hpp
new file mode 100644
index 0000000..98e1f63
--- /dev/null
+++ b/include/nn_search.hpp
@@ -0,0 +1,134 @@
+#pragma once
+
+#include <nanobind/eigen/dense.h>
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/pair.h>
+
+#include <Eigen/Dense>
+#include <iostream>
+#include <nanoflann.hpp>
+#include <taskflow/algorithm/for_each.hpp>
+
+#include "pca.hpp"
+namespace nb = nanobind;
+
+namespace pgeof
+{
+
+/**
+ * Given two point clouds, compute for each point present in one of the point cloud
+ * the N closest points in the other point cloud
+ *
+ * It should be faster than scipy.spatial.KDTree for this task.
+ *
+ * @param data the reference point cloud.
+ * @param query the point cloud used for the queries.
+ * @param knn the number of neighbors to take into account for each point.
+ * @return a pair of nd::array, both of size (n_points x knn), the first one contains the indices of each neighbor, the
+ * second one the square distances between the query point and each of its neighbors.
+ */
+template <typename real_t>
+static std::pair<nb::ndarray<nb::numpy, uint32_t, nb::ndim<2>>, nb::ndarray<nb::numpy, real_t, nb::ndim<2>>>
+    nanoflann_knn_search(RefCloud<real_t> data, RefCloud<real_t> query, const uint32_t knn)
+{
+    using kd_tree_t = nanoflann::KDTreeEigenMatrixAdaptor<RefCloud<real_t>, 3, nanoflann::metric_L2>;
+
+    if (knn > data.rows()) { throw std::invalid_argument("knn size is greater than the data point cloud size"); }
+
+    kd_tree_t          kd_tree(3, data, 10);
+    const Eigen::Index n_points = query.rows();
+    uint32_t*          indices  = new uint32_t[knn * n_points];
+    nb::capsule        owner_indices(indices, [](void* p) noexcept { delete[] (uint32_t*)p; });
+
+    real_t*     sqr_dist = new real_t[knn * n_points];
+    nb::capsule owner_dist(sqr_dist, [](void* p) noexcept { delete[] (real_t*)p; });
+
+    tf::Executor executor;
+    tf::Taskflow taskflow;
+    taskflow.for_each_index(
+        Eigen::Index(0), n_points, Eigen::Index(1),
+        [&](Eigen::Index point_id)
+        {
+            nanoflann::KNNResultSet<real_t, uint32_t, uint32_t> result_set(knn);
+
+            const size_t id = point_id * knn;
+            result_set.init(&indices[id], &sqr_dist[id]);
+            kd_tree.index_->findNeighbors(result_set, query.row(point_id).data());
+        }),
+        tf::StaticPartitioner(0);
+
+    executor.run(taskflow).get();
+
+    const size_t shape[2] = {static_cast<size_t>(n_points), static_cast<size_t>(knn)};
+    return {
+        nb::ndarray<nb::numpy, uint32_t, nb::ndim<2>>(indices, 2, shape, owner_indices),
+        nb::ndarray<nb::numpy, real_t, nb::ndim<2>>(sqr_dist, 2, shape, owner_dist)};
+};
+
+/**
+ * Search for the points within a specified sphere in a point cloud.
+ *
+ * It could be a fallback replacement for FRNN into SuperPointTransformer code base.
+ * It should be faster than scipy.spatial.KDTree for this task.
+ *
+ * @param data the reference point cloud.
+ * @param query the point cloud used for the queries (sphere centers)
+ * @param search_radius the search radius.
+ * @param max_knn the maximum number of neighbors to fetch inside the radius. (Fixing a
+ * reasonable max number of neighbors prevents running OOM for large radius/dense point clouds.)
+ * @return a pair of nd::array, both of size (n_points x knn), the first one contains the 'indices' of each neighbor,
+ * the second one the 'square_distances' between the query point and each neighbor. Point having a number of neighbors <
+ * 'max_knn' inside the 'search_radius' will have their 'indices' and and 'square_distances' filled respectively with
+ * '-1' and 'O' for any missing neighbor.
+ */
+template <typename real_t>
+static std::pair<nb::ndarray<nb::numpy, int32_t, nb::ndim<2>>, nb::ndarray<nb::numpy, real_t, nb::ndim<2>>>
+    nanoflann_radius_search(
+        RefCloud<real_t> data, RefCloud<real_t> query, const real_t search_radius, const uint32_t max_knn)
+{
+    using kd_tree_t = nanoflann::KDTreeEigenMatrixAdaptor<RefCloud<real_t>, 3, nanoflann::metric_L2>;
+
+    if (max_knn > data.rows())
+    {
+        throw std::invalid_argument("max knn size is greater than the data point cloud size");
+    }
+
+    kd_tree_t    kd_tree(3, data, 10);
+    const real_t sq_search_radius = search_radius * search_radius;
+
+    const Eigen::Index n_points = query.rows();
+
+    int32_t*    indices = new int32_t[max_knn * n_points];
+    nb::capsule owner_indices(indices, [](void* p) noexcept { delete[] (int32_t*)p; });
+    std::fill(indices, indices + (max_knn * n_points), -1);
+
+    real_t*     sqr_dist = new real_t[max_knn * n_points];
+    nb::capsule owner_dist(sqr_dist, [](void* p) noexcept { delete[] (real_t*)p; });
+    std::fill(sqr_dist, sqr_dist + (max_knn * n_points), real_t(0.0));
+
+    tf::Executor executor;
+    tf::Taskflow taskflow;
+
+    taskflow.for_each_index(
+        Eigen::Index(0), n_points, Eigen::Index(1),
+        [&](Eigen::Index point_id)
+        {
+            nanoflann::RKNNResultSet<real_t, int32_t, uint32_t> result_set(max_knn, sq_search_radius);
+
+            const size_t id = point_id * max_knn;
+
+            result_set.init(&indices[id], &sqr_dist[id]);
+            kd_tree.index_->findNeighbors(result_set, query.row(point_id).data());
+        },
+        tf::StaticPartitioner(0));
+
+    executor.run(taskflow).get();
+
+    const size_t shape[2] = {static_cast<size_t>(n_points), static_cast<size_t>(max_knn)};
+    return {
+        nb::ndarray<nb::numpy, int32_t, nb::ndim<2>>(indices, 2, shape, owner_indices),
+        nb::ndarray<nb::numpy, real_t, nb::ndim<2>>(sqr_dist, 2, shape, owner_dist)};
+};
+
+}  // namespace pgeof
diff --git a/include/pca.hpp b/include/pca.hpp
new file mode 100644
index 0000000..2bbc382
--- /dev/null
+++ b/include/pca.hpp
@@ -0,0 +1,297 @@
+#pragma once
+
+#include <nanobind/eigen/dense.h>
+
+#include <Eigen/Dense>
+#include <Eigen/Eigenvalues>
+
+namespace nb = nanobind;
+
+namespace pgeof
+{
+
+// Type definitions
+template <typename real_t>
+using PointCloud = Eigen::Matrix<real_t, Eigen::Dynamic, 3, Eigen::RowMajor>;
+
+template <typename real_t>
+using RefCloud = Eigen::Ref<const PointCloud<real_t>>;
+
+template <typename real_t>
+using Vec3 = Eigen::RowVector<real_t, 3>;
+
+template <typename real_t>
+using MatrixCloud = Eigen::Matrix<real_t, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+template <typename real_t>
+using DRefMatrixCloud = nb::DRef<const MatrixCloud<real_t>>;
+
+// epsilon definition, for now same for float an double
+// the eps is meant to stabilize the division when the cloud's 3rd eigenvalue is near 0
+template <typename real_t> constexpr real_t epsilon;
+template <> constexpr float  epsilon<float>  = 1e-3f;
+template <> constexpr double epsilon<double> = 1e-3;
+
+template <typename real_t>
+struct PCAResult
+{
+    Vec3<real_t> val;
+    Vec3<real_t> v0;
+    Vec3<real_t> v1;
+    Vec3<real_t> v2;
+};
+
+// enum of features
+typedef enum EFeatureID
+{
+    Linearity = 0,
+    Planarity,
+    Scattering,
+    VerticalityPGEOF,  // Formula is different from the "classical" formula
+    Normal_x,  // Normal as the third eigenvector
+    Normal_y,
+    Normal_z,
+    Length,
+    Surface,
+    Volume,
+    Curvature,
+    K_optimal,
+    Verticality,  // this is the "classical" verticality
+    Eigentropy
+} EFeatureID;
+
+/**
+ * Given A point cloud compute a PCAResult
+ *
+ * @param cloud the point cloud
+ * @returns A PCAResult
+ */
+template <typename real_t>
+static inline PCAResult<real_t> pca_from_pointcloud(const PointCloud<real_t>& cloud)
+{
+    // Compute the (3, 3) covariance matrix
+    const PointCloud<real_t>          centered_cloud = cloud.rowwise() - cloud.colwise().mean();
+    const Eigen::Matrix<real_t, 3, 3> cov = (centered_cloud.adjoint() * centered_cloud) / real_t(cloud.rows());
+
+    // Compute the eigenvalues and eigenvectors of the covariance
+    Eigen::EigenSolver<Eigen::Matrix<real_t, 3, 3>> es(cov);
+
+    // Sort the values and vectors in order of increasing eigenvalue
+    const auto ev = es.eigenvalues().real();
+
+    std::array<Eigen::Index, 3> indices = {0, 1, 2};
+
+    std::sort(
+        std::begin(indices), std::end(indices), [&](Eigen::Index i1, Eigen::Index i2) { return ev(i1) > ev(i2); });
+
+    Vec3<real_t> val = {
+        (std::max(ev(indices[0]), real_t(0.))), (std::max(ev(indices[1]), real_t(0.))),
+        (std::max(ev(indices[2]), real_t(0.)))};
+    Vec3<real_t> v0 = es.eigenvectors().col(indices[0]).real();
+    Vec3<real_t> v1 = es.eigenvectors().col(indices[1]).real();
+    Vec3<real_t> v2 = es.eigenvectors().col(indices[2]).real();
+
+    // To standardize the orientation of eigenvectors, we choose to enforce all eigenvectors
+    // to be expressed in the Z+ half-space.
+    // Only the third eigenvector (v2) needs to be reoriented because it is the
+    // only one used in further computations.
+    // TODO: In case we want to orient normal, this should be improved
+    if (v2(2) < real_t(0.)) { v2 = real_t(-1.) * v2; }
+    return {val, v0, v1, v2};
+};
+
+/**
+ * Given A point cloud and a CSR definition of the neighboring information for each point, compute a PCAResult
+ *
+ * @param xyz the point cloud
+ * @param nn Integer 1D array. Flattened neighbor indices. Make sure those are all positive,
+ *  '-1' indices will either crash or silently compute incorrect
+ * @param nn_ptr: [n_points+1] Integer 1D array. Pointers wrt 'nn'. More specifically, the neighbors of point 'i'
+ *  are 'nn[nn_ptr[i]:nn_ptr[i + 1]]'
+ * @param i_point the index of the 'central point' or  point
+ * @param k_nnn the number of neighbors to take into account to compute the PCA. It's the caller responsibility
+ * to ensure k_nn won't overflow nn_ptr array.
+ * @returns A PCAResult
+ */
+template <typename real_t, typename index_t>
+static PCAResult<real_t> pca_from_neighborhood(
+    RefCloud<real_t> xyz, const index_t* nn, const index_t* nn_ptr, const size_t i_point, const size_t k_nn)
+{
+    // Initialize the cloud (n_neighbors, 3) matrix holding the
+    // points' neighbors XYZ coordinates
+    PointCloud<real_t> cloud(k_nn, 3);
+    // Recover the neighbors' XYZ coordinates using nn and xyz
+    for (size_t i_nei = 0; i_nei < k_nn; i_nei++)
+    {
+        // Recover the neighbor's position in the xyz vector
+        const Eigen::Index idx_nei = static_cast<Eigen::Index>(nn[nn_ptr[i_point] + i_nei]);
+        // Recover the corresponding xyz coordinates
+        cloud.row(i_nei) = xyz.row(idx_nei);
+    }
+    return pca_from_pointcloud(cloud);
+};
+
+/**
+ * Given a PCA result compute the eigentropy
+ *
+ * This Eigentropy is used for neighborhood size selection in some function
+ * and could be used a an individual feature as well
+ *
+ * @param pca PCAResult
+ * @return the eigentropy
+ */
+template <typename real_t>
+static inline real_t compute_eigentropy(const PCAResult<real_t>& pca)
+{
+    // Compute the eigentropy as defined in:
+    // http://lareg.ensg.eu/labos/matis/pdf/articles_revues/2015/isprs_wjhm_15.pdf
+    const real_t       val_sum = pca.val.sum() + epsilon<real_t>;
+    const Vec3<real_t> e       = pca.val / val_sum;
+    return (-e(0) * std::log(e(0) + epsilon<real_t>) - e(1) * std::log(e(1) + epsilon<real_t>) - e(2) * std::log(e(2) + epsilon<real_t>));
+};
+
+/**
+ * Given a PCA result compute a full set of feature (the initial set of feature r)
+ *
+ * This function intends on mimicking the behavior of original PGEOF feature computation
+ *
+ * @param[in] pca PCAResult
+ * @param[out] feature_results the array of resulting features.
+ */
+template <typename real_t>
+static void compute_features(const PCAResult<real_t>& pca, real_t* features)
+{
+    constexpr real_t sq_eps    = real_t(1e-6);
+    constexpr real_t cub_eps   = real_t(1e-9);
+    constexpr real_t one_third = real_t(1.) / real_t(3.);
+
+    // Compute the dimensionality features. The eps term is meant
+    // to stabilize the division when the cloud's 3rd eigenvalue is
+    // near 0 (points lie in 1D or 2D). Note we take the sqrt of the
+    // eigenvalues since the PCA eigenvalues are homogeneous to m²
+    const real_t val0      = std::sqrt(pca.val(0));
+    const real_t val1      = std::sqrt(pca.val(1));
+    const real_t val2      = std::sqrt(pca.val(2));
+    const real_t val0_fact = real_t(1.0) / (val0 + epsilon<real_t>);
+
+    features[EFeatureID::Normal_x]   = pca.v2(0);
+    features[EFeatureID::Normal_y]   = pca.v2(1);
+    features[EFeatureID::Normal_z]   = pca.v2(2);
+    features[EFeatureID::Linearity]  = (val0 - val1) * val0_fact;
+    features[EFeatureID::Planarity]  = (val1 - val2) * val0_fact;
+    features[EFeatureID::Scattering] = val2 * val0_fact;
+    features[EFeatureID::Length]     = val0;
+    features[EFeatureID::Surface]    = std::sqrt(val0 * val1 + sq_eps);
+    features[EFeatureID::Volume]     = std::pow(val0 * val1 * val2 + cub_eps, one_third);
+    features[EFeatureID::Curvature]  = val2 / (val0 + val1 + val2 + epsilon<real_t>);
+
+    // Compute the verticality. NB we account for the edge case
+    // where all features are 0
+    if (val0 > real_t(0.))
+    {
+        const Vec3<real_t> unary_vector = {
+            pca.val(0) * std::abs(pca.v0(0)) + pca.val(1) * std::abs(pca.v1(0)) + pca.val(2) * std::abs(pca.v2(0)),
+            pca.val(0) * std::abs(pca.v0(1)) + pca.val(1) * std::abs(pca.v1(1)) + pca.val(2) * std::abs(pca.v2(1)),
+            // pca.v2 is already absolute value (positive). but we keep the operation for now
+            // since we can come with our own normal orientation or any other external normal orientation routine
+            pca.val(0) * std::abs(pca.v0(2)) + pca.val(1) * std::abs(pca.v1(2)) + pca.val(2) * std::abs(pca.v2(2))};
+
+        features[EFeatureID::VerticalityPGEOF] = unary_vector(2) / unary_vector.norm();
+    }
+};
+
+/**
+ * Given a PCA result compute only a subset of features.
+ *
+ * This function intends on mimicking the behavior of Jakteristics
+ *
+ * @param[in] pca PCAResult
+ * @param[in] selected_feature a vector of the type of features to compute
+ * @param[out] feature_result the array of resulting features. Result are inserted sequentially the order is defined
+ * by the selected_feature array.
+ */
+template <typename real_t>
+void compute_selected_features(
+    const PCAResult<real_t>& pca, const std::vector<EFeatureID>& selected_feature, real_t* feature_results)
+{
+    // Compute the dimensionality features. The 1e-3 term is meant
+    // to stabilize the division when the cloud's 3rd eigenvalue is
+    // near 0 (points lie in 1D or 2D). Note we take the sqrt of the
+    // eigenvalues since the PCA eigenvalues are homogeneous to m²
+    const real_t val0      = std::sqrt(pca.val(0));
+    const real_t val1      = std::sqrt(pca.val(1));
+    const real_t val2      = std::sqrt(pca.val(2));
+    const real_t val0_fact = real_t(1.0) / (val0 + epsilon<real_t>);
+
+    const auto compute_feature = [val0, val1, val2, val0_fact, &pca](
+                                     const EFeatureID feature_id, const size_t output_id, auto* feature_results)
+    {
+        switch (feature_id)
+        {
+            case EFeatureID::Normal_x:
+                feature_results[output_id] = pca.v2(0);
+                break;
+            case EFeatureID::Normal_y:
+                feature_results[output_id] = pca.v2(1);
+                break;
+            case EFeatureID::Normal_z:
+                feature_results[output_id] = pca.v2(2);
+                break;
+            case EFeatureID::Linearity:
+                feature_results[output_id] = (val0 - val1) * val0_fact;
+                break;
+            case EFeatureID::Planarity:
+                feature_results[output_id] = (val1 - val2) * val0_fact;
+                break;
+            case EFeatureID::Scattering:
+                feature_results[output_id] = val2 * val0_fact;
+                break;
+            case EFeatureID::Length:
+                feature_results[output_id] = val0;
+                break;
+            case EFeatureID::Surface:
+                feature_results[output_id] = std::sqrt(val0 * val1 + 1e-6f);
+                break;
+            case EFeatureID::Volume:
+                feature_results[output_id] = std::pow(
+                    val0 * val1 * val2 + real_t(1e-8),
+                    real_t(1.) / real_t(3.));  // 1e-9 eps is a too small value for float32 so we fallback to 1e-8
+                break;
+            case EFeatureID::Curvature:
+                feature_results[output_id] = val2 / (val0 + val1 + val2 + epsilon<real_t>);
+                break;
+            case EFeatureID::VerticalityPGEOF:
+                // the verticality as defined in PGEOF
+                if (val0 > real_t(0.))
+                {
+                    const Vec3<real_t> unary_vector = {
+                        pca.val(0) * std::abs(pca.v0(0)) + pca.val(1) * std::abs(pca.v1(0)) +
+                            pca.val(2) * std::abs(pca.v2(0)),
+                        pca.val(0) * std::abs(pca.v0(1)) + pca.val(1) * std::abs(pca.v1(1)) +
+                            pca.val(2) * std::abs(pca.v2(1)),
+                        // pca.v2 is already absolute value (positive). But we keep the operation for now
+                        // since we can come with our own normal orientation or any other external normal orientation
+                        // routine
+                        pca.val(0) * std::abs(pca.v0(2)) + pca.val(1) * std::abs(pca.v1(2)) +
+                            pca.val(2) * std::abs(pca.v2(2))};
+
+                    feature_results[output_id] = unary_vector(2) / unary_vector.norm();
+                    // TODO: Jakteristics compute this as feature_results[output_id] = real_t(1.0) -
+                    // std::abs(pca.v2(2));
+                    // It seems to be the most common formula for the verticality in the literature
+                }
+                break;
+            case EFeatureID::Verticality:
+                // The verticality as defined in most of the papers
+                // http://lareg.ensg.eu/labos/matis/pdf/articles_revues/2015/isprs_wjhm_15.pdf
+                feature_results[output_id] = real_t(1.0) - std::abs(pca.v2(2));
+            case EFeatureID::Eigentropy:
+                feature_results[output_id] = compute_eigentropy(pca);
+            default:
+                break;
+        }
+    };
+
+    for (size_t i = 0; i < selected_feature.size(); ++i) { compute_feature(selected_feature[i], i, feature_results); }
+}
+
+}  // namespace pgeof
diff --git a/include/pgeof.hpp b/include/pgeof.hpp
index 5340bc1..2f5203f 100644
--- a/include/pgeof.hpp
+++ b/include/pgeof.hpp
@@ -1,274 +1,369 @@
-#include <cstdio>
+#pragma once
+
+#include <nanobind/eigen/dense.h>
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+
+#include <algorithm>
 #include <cmath>
+#include <cstdio>
 #include <iostream>
+#include <limits>
+#include <taskflow/algorithm/for_each.hpp>
+#include <taskflow/taskflow.hpp>
 #include <vector>
 
-#ifdef PGEOF_WINDOWS
-    #include <ppl.h>
-#endif
+#include "pca.hpp"
 
-#include <eigen3/Eigen/Dense>
-#include <eigen3/Eigen/Eigenvalues>
+namespace nb = nanobind;
 
-#pragma once
-
-namespace ei = Eigen;
+namespace pgeof
+{
 
-struct PCAOutput {
-    std::array<float, 3> val;
-    std::array<float, 3> v0;
-    std::array<float, 3> v1;
-    std::array<float, 3> v2;
-    float eigenentropy;
+namespace log
+{
+/**
+ * Log a progress into the std::cout. It has multiple quirks (std::cout usage, race condition)
+ * but it is mostly used for debugging purposes.
+ *
+ * @param progress_count a reference to the current count. it is incremented by one at each call of this function
+ * @param progress_total the total expected count of iteration
+ */
+static inline void progress(size_t& progress_count, const size_t progress_total)
+{
+    ++progress_count;
+    // Print progress
+    // NB: when in parallel progress_count behavior is undefined, but
+    // gives a good indication of progress
+    if (progress_count % 10000 == 0)
+    {
+        std::cout << progress_count << "% done          \r" << std::flush;
+        std::cout << std::ceil(progress_count * 100 / progress_total) << "% done          \r" << std::flush;
+    }
 };
-typedef ei::Matrix<float, 3, 3> Matrix3f;
-typedef ei::Matrix<float, 3, 1> Vector3f;
-typedef ei::Matrix<float, ei::Dynamic, 3> PointCloud;
 
-PCAOutput neighborhood_pca(
-    float *xyz, uint32_t *nn, const uint32_t *nn_ptr, std::size_t i_point,
-    std::size_t k_nn)
+/**
+ * flush the logger
+ */
+static inline void flush() { std::cout << std::endl; };
+}  // namespace log
+
+/**
+ * Compute a set of geometric features for a point cloud from a precomputed list of neighbors.
+ *
+ *  * The following features are computed:
+ * - linearity
+ * - planarity
+ * - scattering
+ * - verticality
+ * - normal vector (oriented towards positive z-coordinates)
+ * - length
+ * - surface
+ * - volume
+ * - curvature
+ *
+ * @param xyz The point cloud.
+ * @param nn Integer 1D array. Flattened neighbor indices. Make sure those are all positive,
+ * '-1' indices will either crash or silently compute incorrect features.
+ * @param nn_ptr: [n_points+1] Integer 1D array. Pointers wrt 'nn'. More specifically, the neighbors of point 'i'
+ * are 'nn[nn_ptr[i]:nn_ptr[i + 1]]'.
+ * @param k_min Minimum number of neighbors to consider for features computation. If a point has less,
+ * it features will be a set of '0' values.
+ * @param verbose Whether computation progress should be printed out
+ * @return the geometric features associated with each point's neighborhood in a (num_points, features_count) nd::array.
+ */
+template <typename real_t = float, const size_t feature_count = 11>
+static nb::ndarray<nb::numpy, real_t, nb::shape<nb::any, feature_count>> compute_geometric_features(
+    RefCloud<real_t> xyz, nb::ndarray<const uint32_t, nb::ndim<1>> nn, nb::ndarray<const uint32_t, nb::ndim<1>> nn_ptr,
+    const size_t k_min, const bool verbose)
 {
-    // Initialize the cloud (n_neighbors + 1, 3) matrix holding the
-    // points' neighbors XYZ coordinates
-    PointCloud cloud(k_nn, 3);
+    if (k_min < 1) { throw std::invalid_argument("k_min should be > 1"); }
+    // Each point can be treated in parallel
+    const size_t    n_points    = nn_ptr.size() - 1;  // number of points is not determined by xyz
+    size_t          s_point     = 0;
+    const uint32_t* nn_data     = nn.data();
+    const uint32_t* nn_ptr_data = nn_ptr.data();
 
-    // Recover the neighbors' XYZ coordinates using nn and xyz
-    std::size_t idx_nei;
-    for (std::size_t i_nei = 0; i_nei < k_nn; i_nei++)
-    {
-        // Recover the neighbor's position in the xyz vector
-        idx_nei = nn[nn_ptr[i_point] + i_nei];
+    real_t*     features = (real_t*)calloc(n_points * feature_count, sizeof(real_t));
+    nb::capsule owner_features(features, [](void* f) noexcept { delete[] (real_t*)f; });
 
-        // Recover the corresponding xyz coordinates
-        cloud(i_nei, 0) = xyz[3 * idx_nei];
-        cloud(i_nei, 1) = xyz[3 * idx_nei + 1];
-        cloud(i_nei, 2) = xyz[3 * idx_nei + 2];
-    }
+    tf::Executor executor;
+    tf::Taskflow taskflow;
+    taskflow.for_each_index(
+        size_t(0), size_t(n_points), size_t(1),
+        [&](size_t i_point)
+        {
+            if (verbose) log::progress(s_point, n_points);
 
-    // Compute the (3, 3) covariance matrix
-    ei::MatrixXf centered_cloud = cloud.rowwise() - cloud.colwise().mean();
-    ei::Matrix3f cov =
-        (centered_cloud.adjoint() * centered_cloud) / float(k_nn);
-
-    // Compute the eigenvalues and eigenvectors of the covariance
-    ei::EigenSolver<Matrix3f> es(cov);
-
-    // Sort the values and vectors in order of increasing eigenvalue
-    std::vector<float> ev = {
-        es.eigenvalues()[0].real(),
-        es.eigenvalues()[1].real(),
-        es.eigenvalues()[2].real()};
-    std::array<std::size_t, 3> indices = {0, 1, 2};
-    std::sort(
-        std::begin(indices),
-        std::end(indices),
-        [&](std::size_t i1, std::size_t i2) { return ev[i1] > ev[i2]; } );
-    std::array<float,3> val = {
-        (std::max(ev[indices[0]],0.f)),
-        (std::max(ev[indices[1]],0.f)),
-        (std::max(ev[indices[2]],0.f))};
-    std::array<float,3> v0 = {
-        es.eigenvectors().col(indices[0])(0).real(),
-        es.eigenvectors().col(indices[0])(1).real(),
-        es.eigenvectors().col(indices[0])(2).real()};
-    std::array<float,3>  v1 = {
-        es.eigenvectors().col(indices[1])(0).real(),
-        es.eigenvectors().col(indices[1])(1).real(),
-        es.eigenvectors().col(indices[1])(2).real()};
-    std::array<float,3>  v2 = {
-        es.eigenvectors().col(indices[2])(0).real(),
-        es.eigenvectors().col(indices[2])(1).real(),
-        es.eigenvectors().col(indices[2])(2).real()};
-
-    // To standardize the orientation of eigenvectors, we choose to
-    // enforce all eigenvectors to be expressed in the Z+ half-space
-    if (v0[2] < 0)
-    {
-        v0[0] = -v0[0];
-        v0[1] = -v0[1];
-        v0[2] = -v0[2];
-    }
+            // Recover the points' total number of neighbors
+            const size_t k_nn = static_cast<size_t>(nn_ptr_data[i_point + 1] - nn_ptr_data[i_point]);
+
+            // If the cloud has less than k_min point, continue
+            if (k_nn >= k_min)
+            {
+                const PCAResult<real_t> pca = pca_from_neighborhood(xyz, nn_data, nn_ptr_data, i_point, k_nn);
+                compute_features(pca, &features[i_point * feature_count]);
+            }
+        },
+        tf::StaticPartitioner(0));
+    executor.run(taskflow).get();
 
-    if (v1[2] < 0)
+    // Final print to start on a new line
+    if (verbose) log::flush();
+    const size_t shape[2] = {n_points, feature_count};
+    return nb::ndarray<nb::numpy, real_t, nb::shape<nb::any, feature_count>>(features, 2, shape, owner_features);
+}
+/**
+ * Convenience function that check that scales are well ordered in increasing order.
+ *
+ * @param k_scales the list of scale size (number of neighbors).
+ */
+static bool check_scales(const std::vector<uint32_t>& k_scales)
+{
+    uint32_t previous_scale = 1;  // minimal admissible k_min value is 1
+    for (const auto& current_scale : k_scales)
     {
-        v1[0] = -v1[0];
-        v1[1] = -v1[1];
-        v1[2] = -v1[2];
+        if (current_scale < previous_scale) { return false; }
+        previous_scale = current_scale;
     }
+    return true;
+}
 
-    if (v2[2] < 0)
+/**
+ * Compute a set of geometric features for a point cloud in a multiscale fashion.
+ *
+ * The following features are computed:
+ * - linearity
+ * - planarity
+ * - scattering
+ * - verticality
+ * - normal vector (oriented towards positive z-coordinates)
+ * - length
+ * - surface
+ * - volume
+ * - curvature
+ *
+ * @param xyz The point cloud
+ * @param nn Integer 1D array. Flattened neighbor indices. Make sure those are all positive,
+ *  '-1' indices will either crash or silently compute incorrect features.
+ * @param nn_ptr: [n_points+1] Integer 1D array. Pointers wrt 'nn'. More specifically, the neighbors of point 'i'
+ *  are 'nn[nn_ptr[i]:nn_ptr[i + 1]]'.
+ * @param k_scale Array of number of neighbors to consider for features computation. If a at a given scale, a point has
+ * less features will be a set of '0' values.
+ * @param verbose Whether computation progress should be printed out
+ * @return Geometric features associated with each point's neighborhood in a (num_points, features_count, n_scales)
+ * nd::array
+ */
+template <typename real_t, const size_t feature_count = 11>
+static nb::ndarray<nb::numpy, real_t, nb::shape<nb::any, nb::any, feature_count>> compute_geometric_features_multiscale(
+    RefCloud<real_t> xyz, nb::ndarray<const uint32_t, nb::ndim<1>> nn, nb::ndarray<const uint32_t, nb::ndim<1>> nn_ptr,
+    const std::vector<uint32_t>& k_scales, const bool verbose)
+{
+    if (!check_scales(k_scales))
     {
-        v2[0] = -v2[0];
-        v2[1] = -v2[1];
-        v2[2] = -v2[2];
+        throw std::invalid_argument("k_scales should be > 1 and sorted in ascending order");
     }
+    const size_t    n_points    = nn_ptr.size() - 1;  // number of points is not determined by xyz
+    const size_t    n_scales    = k_scales.size();
+    size_t          s_point     = 0;
+    const uint32_t* nn_data     = nn.data();
+    const uint32_t* nn_ptr_data = nn_ptr.data();
 
-    // Compute the eigenentropy as defined in:
-    // http://lareg.ensg.eu/labos/matis/pdf/articles_revues/2015/isprs_wjhm_15.pdf
-    const float epsilon = 0.001f;
-    const float val_sum = val[0] + val[1] + val[2] + epsilon;
-    const float e0 = val[0] / val_sum;
-    const float e1 = val[1] / val_sum;
-    const float e2 = val[2] / val_sum;
-    const float eigenentropy = - e0 * std::log(e0 + epsilon) - e1 * std::log(e1 + epsilon)
-        - e2 * std::log(e2 + epsilon);
-
-    return {val, v0, v1, v2, eigenentropy};
-}
+    real_t*     features = (real_t*)calloc(n_points * n_scales * feature_count, sizeof(real_t));
+    nb::capsule owner_features(features, [](void* f) noexcept { delete[] (real_t*)f; });
+
+    // Each point can be treated in parallel
+    tf::Executor executor;
+    tf::Taskflow taskflow;
+    taskflow.for_each_index(
+        size_t(0), size_t(n_points), size_t(1),
+        [&](size_t i_point)
+        {
+            if (verbose) log::progress(s_point, n_points);
+            // Recover the points' total number of neighbors
+            const size_t k_nn = static_cast<size_t>(nn_ptr_data[i_point + 1] - nn_ptr_data[i_point]);
+
+            for (size_t i_scale = 0; i_scale < n_scales; ++i_scale)
+            {
+                const size_t knn_scale = static_cast<size_t>(k_scales[i_scale]);
+
+                if (k_nn < knn_scale)
+                    break;  // we assume scales are stored in increasing order,
+                            // so we could do an early break in case of k_nn <
+                            // knn_scale
+                const PCAResult<real_t> pca = pca_from_neighborhood(xyz, nn_data, nn_ptr_data, i_point, knn_scale);
+                compute_features(pca, &features[(i_point * n_scales + i_scale) * feature_count]);
+            }
+        },
+        tf::StaticPartitioner(0));
+
+    executor.run(taskflow).get();
 
+    // Final print to start on a new line
+    if (verbose) log::flush();
+
+    const size_t shape[3] = {n_points, n_scales, feature_count};
+    return nb::ndarray<nb::numpy, real_t, nb::shape<nb::any, nb::any, feature_count>>(
+        features, 3, shape, owner_features);
+}
 
-void compute_geometric_features(
-    float *xyz, uint32_t *nn, const uint32_t *nn_ptr, int n_points,
-    float *features, int k_min, int k_step, int k_min_search, bool verbose)
+/**
+ * Compute a set of geometric features for a point cloud using the optimal neighborhood selection described in
+ * http://lareg.ensg.eu/labos/matis/pdf/articles_revues/2015/isprs_wjhm_15.pdf
+ *
+ *  * The following features are computed:
+ * - linearity
+ * - planarity
+ * - scattering
+ * - verticality
+ * - normal vector (oriented towards positive z-coordinates)
+ * - length
+ * - surface
+ * - volume
+ * - curvature
+ * - optimal_nn
+ *
+ * @param xyz The point cloud
+ * @param nn Integer 1D array. Flattened neighbor indices. Make sure those are all positive,
+ *  '-1' indices will either crash or silently compute incorrect features.
+ * @param nn_ptr: [n_points+1] Integer 1D array. Pointers wrt 'nn'. More specifically, the neighbors of point 'i'
+ *  are 'nn[nn_ptr[i]:nn_ptr[i + 1]]'.
+ * @param k_min Minimum number of neighbors to consider for features computation. If a point has less,
+ * its features will be a set of '0' values.
+ * @param k_step Step size to take when searching for the optimal neighborhood, size for each point following
+ * Weinmann, 2015
+ * @param k_min_search Minimum neighborhood size at which to start when searching for the optimal neighborhood size for
+ each point. It is advised to use a value of 10 or higher, for geometric features robustness.
+ * @param verbose Whether computation progress should be printed out
+ * @return Geometric features associated with each point's neighborhood in a (num_points, features_count) nd::array
+ */
+template <typename real_t, const size_t feature_count = 12>
+static nb::ndarray<nb::numpy, real_t, nb::shape<nb::any, feature_count>> compute_geometric_features_optimal(
+    RefCloud<real_t> xyz, nb::ndarray<const uint32_t, nb::ndim<1>> nn, nb::ndarray<const uint32_t, nb::ndim<1>> nn_ptr,
+    const uint32_t k_min, const uint32_t k_step, const uint32_t k_min_search, const bool verbose)
 {
+    if (k_min < 1 && k_min_search < 1) { throw std::invalid_argument("k_min and k_min_search should be > 1"); }
     // Each point can be treated in parallel
-    std::size_t s_point = 0;
-
-    #ifdef PGEOF_WINDOWS
-    concurrency::parallel_for(std::size_t(0), std::size_t(n_points), [&](std::size_t i_point)
-    #else
-    #pragma omp parallel for schedule(static)
-    for (std::size_t i_point = 0; i_point < n_points; i_point++)
-    #endif
-    {
-        // Recover the points' total number of neighbors
-        const std::size_t k_nn = nn_ptr[i_point + 1] - nn_ptr[i_point];
+    const size_t    n_points    = nn_ptr.size() - 1;  // number of points is not determined by xyz
+    size_t          s_point     = 0;
+    const uint32_t* nn_data     = nn.data();
+    const uint32_t* nn_ptr_data = nn_ptr.data();
 
-        // If the cloud has only one point, populate the final feature
-        // vector with zeros and continue
-        if (k_nn < k_min ||  k_nn <= 0)
-        {
-            features[i_point * 12 + 0]  = 0.f;
-            features[i_point * 12 + 1]  = 0.f;
-            features[i_point * 12 + 2]  = 0.f;
-            features[i_point * 12 + 3]  = 0.f;
-            features[i_point * 12 + 4]  = 0.f;
-            features[i_point * 12 + 5]  = 0.f;
-            features[i_point * 12 + 6]  = 0.f;
-            features[i_point * 12 + 7]  = 0.f;
-            features[i_point * 12 + 8]  = 0.f;
-            features[i_point * 12 + 9]  = 0.f;
-            features[i_point * 12 + 10] = 0.f;
-            features[i_point * 12 + 11] = 0.f;
-            s_point++;
-            #ifdef PGEOF_WINDOWS
-            return;
-            #else
-            continue;
-            #endif
-        }
-
-        // Compute the PCA for neighborhoods of increasing sizes. The
-        // neighborhood size with the lowest eigenentropy will be kept.
-        // Do not search optimal neighborhood size if k_step < 1
-        PCAOutput pca;
-        std::size_t k_optimal =  k_nn;
-        if (k_step < 1)
-        {
-            pca = neighborhood_pca(xyz, nn, nn_ptr, i_point, k_nn);
-        }
-        else
+    real_t*     features = (real_t*)calloc(n_points * feature_count, sizeof(real_t));
+    nb::capsule owner_features(features, [](void* f) noexcept { delete[] (real_t*)f; });
+
+    tf::Executor executor;
+    tf::Taskflow taskflow;
+    taskflow.for_each_index(
+        size_t(0), size_t(n_points), size_t(1),
+        [&](size_t i_point)
         {
-            // We do not want to search too-small neighborhoods, since
-            // their eigenentropy is likely to be small despite having
-            // unreliable or noisy geometric features. Hence, we start
-            // searching at k_min_search
-            std::size_t k0 = std::min(std::max(k_min, k_min_search), int(k_nn));
+            if (verbose) log::progress(s_point, n_points);
 
-            for (std::size_t k = k0; k < k_nn + 1; ++k)
-            {
-                // Only evaluate the neighborhood's PCA every 'k_step'
-                // and at the boundary values: k0 and k_nn
-                if ((k > k0) && (k % k_step != 0) && (k != k_nn))
-                {
-                    continue;
-                }
+            // Recover the points' total number of neighbors
+            const size_t k_nn = static_cast<size_t>(nn_ptr_data[i_point + 1] - nn_ptr_data[i_point]);
 
-                // Actual PCA computation on the k-neighborhood
-                const PCAOutput pca_k = neighborhood_pca(xyz, nn, nn_ptr, i_point, k);
+            // Process only if the cloud has the required number of point
+            if (k_nn >= k_min && k_nn >= k_min_search)
+            {
+                size_t k0 = std::min(std::max(static_cast<size_t>(k_min), static_cast<size_t>(k_min_search)), k_nn);
 
-                // Keep track of the optimal neighborhood size with the
-                // lowest eigenentropy
-                if ((k == k0) || (pca_k.eigenentropy < pca.eigenentropy))
+                PCAResult<real_t> pca_optimal;
+                real_t            eigenentropy_optimal = real_t(1.0);
+                size_t            k_optimal            = k_nn;
+                for (size_t k = k0; k <= k_nn; ++k)
                 {
-                    pca = pca_k;
-                    k_optimal = k;
-                    continue;
+                    // Only evaluate the neighborhood's PCA every 'k_step'
+                    // and at the boundary values: k0 and k_nn
+                    if ((k > k0) && (k % k_step != 0) && (k != k_nn)) { continue; }
+
+                    const PCAResult<real_t> pca          = pca_from_neighborhood(xyz, nn_data, nn_ptr_data, i_point, k);
+                    const real_t            eigenentropy = compute_eigentropy(pca);
+                    // Keep track of the optimal neighborhood size with the
+                    // lowest eigenentropy
+                    if ((k == k0) || (eigenentropy < eigenentropy_optimal))
+                    {
+                        eigenentropy_optimal = eigenentropy;
+                        k_optimal            = k;
+                        pca_optimal          = pca;
+                    }
                 }
+                compute_features(pca_optimal, &features[i_point * feature_count]);
+                // Add best nn
+                features[i_point * feature_count + 11] = real_t(k_optimal);
             }
-        }
-
-        // Recover the eigenvalues and eigenvectors from the PCA
-        const std::array<float, 3> & val = pca.val;
-        const std::array<float, 3> & v0 = pca.v0;
-        const std::array<float, 3> & v1 = pca.v1;
-        const std::array<float, 3> & v2 = pca.v2;
-
-        // Compute the dimensionality features. The 1e-3 term is meant
-        // to stabilize the division when the cloud's 3rd eigenvalue is
-        // near 0 (points lie in 1D or 2D). Note we take the sqrt of the
-        // eigenvalues since the PCA eigenvaluess are homogeneous to m²
-        const float val0       = std::sqrt(val[0]);
-        const float val1       = std::sqrt(val[1]);
-        const float val2       = std::sqrt(val[2]);
-        const float linearity  = (val0 - val1) / (val0 + 1e-3f);
-        const float planarity  = (val1 - val2) / (val0 + 1e-3f);
-        const float scattering = val2 / (val0 + 1e-3f);
-        const float length     = val0;
-        const float surface    = std::sqrt(val0 * val1 + 1e-6f);
-        const float volume     = std::pow(val0 * val1 * val2 + 1e-9f, 1.f / 3.f);
-        const float curvature  = val2 / (val0 + val1 + val2 + 1e-3f);
-
-        // Compute the verticality. NB we account for the edge case
-        // where all features are 0
-        float verticality = 0.f;
-        if (val0 > 0)
-        {
-            const std::array<float,3> unary_vector = {
-                val[0] * std::abs(v0[0]) + val[1] * std::abs(v1[0]) + val[2] * std::abs(v2[0]),
-                val[0] * std::abs(v0[1]) + val[1] * std::abs(v1[1]) + val[2] * std::abs(v2[1]),
-                val[0] * std::abs(v0[2]) + val[1] * std::abs(v1[2]) + val[2] * std::abs(v2[2])};
-            const float norm = std::sqrt(
-                unary_vector[0] * unary_vector[0]
-                + unary_vector[1] * unary_vector[1]
-                + unary_vector[2] * unary_vector[2]);
-            verticality = unary_vector[2] / norm;
-        }
-
-        // Populate the final feature vector
-        features[i_point * 12 + 0]  = linearity;
-        features[i_point * 12 + 1]  = planarity;
-        features[i_point * 12 + 2]  = scattering;
-        features[i_point * 12 + 3]  = verticality;
-        features[i_point * 12 + 4]  = v2[0];
-        features[i_point * 12 + 5]  = v2[1];
-        features[i_point * 12 + 6]  = v2[2];
-        features[i_point * 12 + 7]  = length;
-        features[i_point * 12 + 8]  = surface;
-        features[i_point * 12 + 9]  = volume;
-        features[i_point * 12 + 10] = curvature;
-        features[i_point * 12 + 11] = float(k_optimal);
-
-        // Print progress
-        // NB: when in parallel s_point behavior is undefined, but gives
-        // a good indication of progress
-        s_point++;
-        if (verbose && s_point % 10000 == 0)
+        },
+        tf::StaticPartitioner(0));
+
+    executor.run(taskflow).get();
+
+    if (verbose) log::flush();
+
+    const size_t shape[2] = {n_points, feature_count};
+    return nb::ndarray<nb::numpy, real_t, nb::shape<nb::any, feature_count>>(features, 2, shape, owner_features);
+}
+
+/**
+ * Compute a selected set of geometric features for a point cloud via radius search.
+ *
+ * This function aims to mimic the behavior of jakteristics and provide an efficient way
+ * to compute a limited set of features.
+ *
+ * @param xyz The point cloud
+ * @param search_radius the search radius.
+ * @param max_knn the maximum number of neighbors to fetch inside the radius. The central point is included. Fixing a
+ * reasonable max number of neighbors prevents running OOM for large radius/dense point clouds.
+ * @param selected_features the list of selected features. See pgeof::EFeatureID
+ * @return Geometric features associated with each point's neighborhood in a (num_points, features_count) nd::array
+ */
+template <typename real_t>
+static nb::ndarray<nb::numpy, real_t, nb::shape<nb::any, nb::any>> compute_geometric_features_selected(
+    RefCloud<real_t> xyz, const real_t search_radius, const uint32_t max_knn,
+    const std::vector<EFeatureID>& selected_features)
+{
+    using kd_tree_t = nanoflann::KDTreeEigenMatrixAdaptor<RefCloud<real_t>, 3, nanoflann::metric_L2>;
+    // TODO: where knn < num of points
+    kd_tree_t    kd_tree(3, xyz, 10);
+    const real_t sq_search_radius = search_radius * search_radius;
+    const size_t feature_count    = selected_features.size();
+
+    const Eigen::Index n_points = xyz.rows();
+    real_t*            features = (real_t*)calloc(n_points * feature_count, sizeof(real_t));
+    std::fill(features, features + (n_points * feature_count), real_t(0.0));
+    nb::capsule owner_features(features, [](void* f) noexcept { delete[] (real_t*)f; });
+
+    tf::Executor executor;
+    tf::Taskflow taskflow;
+    taskflow.for_each_index(
+        Eigen::Index(0), n_points, Eigen::Index(1),
+        [&](Eigen::Index point_id)
         {
-            std::cout << s_point << "% done          \r" << std::flush;
-            std::cout << ceil(s_point * 100 / n_points) << "% done          \r" << std::flush;
-        }
-    }
-    #ifdef PGEOF_WINDOWS
-    );
-    #endif
+            nanoflann::RKNNResultSet<real_t, int32_t, uint32_t> result_set(max_knn, sq_search_radius);
 
-    // Final print to start on a new line
-    if (verbose)
-    {
-        std::cout << std::endl;
-    }
+            std::vector<int32_t> indices(max_knn, -1);  // TODO: maybe int64
+            std::vector<real_t>  sqr_dist(max_knn, real_t(0.0));
+
+            result_set.init(&indices[0], &sqr_dist[0]);
+            kd_tree.index_->findNeighbors(result_set, xyz.row(point_id).data());
+
+            std::vector<int32_t> final_indices;
+            for (size_t i = 0; i < max_knn; ++i)
+            {
+                if (indices[i] == -1) { break; }
+                final_indices.push_back(indices[i]);
+            }
+
+            if (final_indices.size() > 2)
+            {
+                const PointCloud<real_t> cloud = xyz(final_indices, Eigen::all);
+                const PCAResult<real_t>  pca   = pca_from_pointcloud(cloud);
+                compute_selected_features(pca, selected_features, &features[point_id * feature_count]);
+            }
+        },
+        tf::StaticPartitioner(0));
+    executor.run(taskflow).get();
 
-    return;
+    const size_t shape[2] = {static_cast<size_t>(n_points), feature_count};
+    return nb::ndarray<nb::numpy, real_t, nb::shape<nb::any, nb::any>>(features, 2, shape, owner_features);
 }
+}  // namespace pgeof
diff --git a/install.sh b/install.sh
deleted file mode 100755
index f2e75f3..0000000
--- a/install.sh
+++ /dev/null
@@ -1,60 +0,0 @@
-#!/bin/bash
-
-# Local variables
-PROJECT_NAME=pgeof
-PYTHON=3.8
-
-
-# Installation script for Anaconda3 environments
-echo "#############################################"
-echo "#                                           #" 
-echo "#          Point Geometric Features         #"
-echo "#                 Installer                 #"
-echo "#                                           #" 
-echo "#############################################"
-echo
-echo
-
-
-echo "_______________ Prerequisites _______________"
-echo "  - conda"
-echo
-
-
-echo "____________ Pick conda install _____________"
-echo
-# Recover the path to conda on your machine
-CONDA_DIR=`realpath ~/miniconda3`
-if (test -z $CONDA_DIR) || [ ! -d $CONDA_DIR ]
-then
-  CONDA_DIR=`realpath ~/anaconda3`
-fi
-
-while (test -z $CONDA_DIR) || [ ! -d $CONDA_DIR ]
-do
-    echo "Could not find conda at: "$CONDA_DIR
-    read -p "Please provide you conda install directory: " CONDA_DIR
-    CONDA_DIR=`realpath $CONDA_DIR`
-done
-
-echo "Using conda conda found at: ${CONDA_DIR}/etc/profile.d/conda.sh"
-source ${CONDA_DIR}/etc/profile.d/conda.sh
-echo
-echo
-
-
-echo "________________ Installation _______________"
-echo
-
-# Create a conda environment from yml
-conda create --name $PROJECT_NAME python=$PYTHON -y
-
-# Activate the env
-source ${CONDA_DIR}/etc/profile.d/conda.sh  
-conda activate ${PROJECT_NAME}
-
-# Dependencies
-conda install -c anaconda pip 
-conda install -c omnia eigen3 -y
-export EIGEN_LIB_PATH="$CONDA_PREFIX/include"   # Add Eigen includes to the path for the C++ compiler
-python -m pip install .
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index b11ed50..d5527f6 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,10 +1,10 @@
 [build-system]
-requires = ["wheel", "setuptools", "oldest-supported-numpy"]
-build-backend = "setuptools.build_meta"
+requires = ["scikit-build-core >=0.4.3", "nanobind >= 1.9.2"]
+build-backend = "scikit_build_core.build"
 
 [project]
 name = "pgeof"
-version = "0.0.2"
+version = "0.1.0"
 readme = "README.md"
 description = "Compute the geometric features associated with each point's neighborhood:"
 requires-python = ">=3.8,<3.13"
@@ -19,14 +19,56 @@ classifiers = [
   "Programming Language :: Python",
   "Topic :: Scientific/Engineering"
 ]
-
 dependencies = [
-    "numpy"
+  "numpy >= 1.7, < 2"
 ]
 
 [project.urls]
 homepage = "https://github.com/drprojects/point_geometric_features"
 repository = "https://github.com/drprojects/point_geometric_features"
 
+[tool.scikit-build]
+# Protect the configuration against future changes in scikit-build-core
+minimum-version = "0.4"
+
+# Setuptools-style build caching in a local directory
+build-dir = "build/{wheel_tag}"
+
+cmake.build-type = "Release"
+
+[tool.tox]
+legacy_tox_ini = """
+[tox]
+
+[gh-actions]
+python =
+    3.9: py39
+    3.10: py310
+    3.11: py311
+    3.12: py312
+
+[testenv]
+deps =
+  pytest >= 7.4
+  pytest-benchmark ~= 4.0
+  numpy >= 1.7, < 2
+  scipy
+  jakteristics;platform_system=="Windows" or platform_system=="Linux" 
+commands = pytest --basetemp="{envtmpdir}" {posargs}
+
+[testenv:bench]
+# globs/wildcards do not work with tox
+commands = pytest --basetemp="{envtmpdir}" {posargs:tests/bench_knn.py tests/bench_jakteristics.py} 
+
+[testenv:ls]
+skipbuild=true
+commands = ls {envtmpdir}
+"""
+
 [tool.cibuildwheel]
-build = "cp3{8,9,10,11,12}-*"
+build = "cp3{9,10,11,12}-*"
+archs = ["auto64"] # limit to 64bits builds
+
+# Needed for full C++17 support
+[tool.cibuildwheel.macos.environment]
+MACOSX_DEPLOYMENT_TARGET = "10.14"
\ No newline at end of file
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 6e4d4a6..0000000
--- a/setup.py
+++ /dev/null
@@ -1,47 +0,0 @@
-# ----------------------------------------------------------------------#
-#      setuptools setup script for compiling pgeof extensions           #
-# ----------------------------------------------------------------------#
-""" 
-Loic Landrieu 2023
-"""
-
-from setuptools import setup, Extension
-import numpy
-import platform
-import os
-
-# targets and compile options
-name = "pgeof"
-
-include_dirs = [numpy.get_include(), "include"]  # find the Numpy headers
-
-EIGEN_LIB_PATH = os.environ.get("EIGEN_LIB_PATH", None)
-
-if EIGEN_LIB_PATH is not None:
-    include_dirs.append(EIGEN_LIB_PATH)
-
-
-# Compilation and linkage options
-if platform.system() == "Windows":
-    extra_compile_args = ["/DPGEOF_WINDOWS"]
-    extra_link_args = []
-elif platform.system() == "Linux":
-    extra_compile_args = ["-std=c++11", "-fopenmp"]
-    extra_link_args = ["-lgomp"]
-elif platform.system() == "Darwin":
-    extra_compile_args = ["-std=c++11", "-fopenmp"]
-    extra_link_args = ["-lomp"]
-else:
-    raise NotImplementedError("OS not yet supported.")
-
-
-#  Compilation
-mod = Extension(
-    name,
-    # list source files
-    ["./src/pgeof_cpy.cpp"],
-    include_dirs=include_dirs,
-    extra_compile_args=extra_compile_args,
-    extra_link_args=extra_link_args)
-
-setup(name=name, ext_modules=[mod])
diff --git a/src/pgeof/__init__.py b/src/pgeof/__init__.py
new file mode 100644
index 0000000..d351914
--- /dev/null
+++ b/src/pgeof/__init__.py
@@ -0,0 +1,9 @@
+from .pgeof_ext import (
+    EFeatureID,
+    compute_features,
+    compute_features_multiscale,
+    compute_features_optimal,
+    knn_search,
+    radius_search,
+    compute_features_selected
+)
diff --git a/src/pgeof_cpy.cpp b/src/pgeof_cpy.cpp
deleted file mode 100644
index 6ca7cc1..0000000
--- a/src/pgeof_cpy.cpp
+++ /dev/null
@@ -1,159 +0,0 @@
-#include <cstdint>
-#include <cstdio>
-#define PY_SSIZE_T_CLEAN
-#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
-#include <Python.h>
-#include <numpy/arrayobject.h>
-
-#include "pgeof.hpp"
-
-/* template for handling several index types in compute_geometric_features */
-static PyObject* pgeof(
-    PyArrayObject *py_xyz, PyArrayObject * py_nn, PyArrayObject * py_nn_ptr,
-    int k_min, int k_step, int k_min_search, bool verbose)
-{
-    //convert from python to arrays
-    float * c_xyz = (float*) PyArray_DATA(py_xyz);
-    uint32_t * c_nn =  (uint32_t*) PyArray_DATA(py_nn);
-    uint32_t * c_nn_ptr= (uint32_t*) PyArray_DATA(py_nn_ptr);
-    int n_points = PyArray_DIMS(py_nn_ptr)[0] - 1;
-    int num_features = 12;
-
-    //prepare output
-    npy_intp size_of_feature[] = {n_points, num_features};
-    PyArrayObject* py_features = (PyArrayObject*) PyArray_Zeros(
-        2, size_of_feature, PyArray_DescrFromType(NPY_FLOAT32), 0);
-    float *features = (float*) PyArray_DATA(py_features);
-
-    compute_geometric_features(
-        c_xyz, c_nn, c_nn_ptr, n_points, features, k_min, k_step, k_min_search, verbose);
-
-    return Py_BuildValue("N", py_features);
-}
-
-
-/* actual interface*/
-static PyObject* pgeof_cpy(PyObject* self, PyObject *args, PyObject *kwargs)
-{   (void) self; // suppress unused parameter warning
-
-    /* inputs  */
-    PyArrayObject *xyz;
-    PyArrayObject *nn;
-    PyArrayObject *nn_ptr;
-    int k_min = 1;
-    int k_step = -1;
-    int k_min_search = 10;
-    bool verbose = false;
-
-    // Build variable names used for input args + kwargs parsing
-    static char *keywords[] = {
-        (char*)"xyz",
-        (char*)"nn",
-        (char*)"nn_ptr",
-        (char*)"k_min",
-        (char*)"k_step",
-        (char*)"k_min_search",
-        (char*)"verbose",
-        NULL};
-
-    /* parse the input, from Python Object to C PyArray */
-//    if(!PyArg_ParseTuple(args, "OOOiiii", &xyz, &nn, &nn_ptr, &k_min, &k_step, &k_min_search, &verbose)){
-//        return NULL;
-//    }
-    if (!PyArg_ParseTupleAndKeywords(
-        args,
-        kwargs,
-        "OOO|iiip",
-        keywords,
-        &xyz,
-        &nn,
-        &nn_ptr,
-        &k_min,
-        &k_step,
-        &k_min_search,
-        &verbose))
-    {
-        return NULL;
-    }
-
-    PyObject* PyReturn = pgeof(xyz, nn, nn_ptr, k_min, k_step, k_min_search, verbose);
-    return PyReturn;
-    }
-
-static const char* pgeof_doc =
-    "Compute the geometric features associated with each point's\n"
-    "neighborhood. The following features are computed:\n"
-    "  - linearity\n"
-    "  - planarity\n"
-    "  - scattering\n"
-    "  - verticality\n"
-    "  - normal vector (oriented towards positive z-coordinates)\n"
-    "  - length\n"
-    "  - surface\n"
-    "  - volume\n"
-    "  - curvature\n"
-    "  - optimal neighborhood size\n\n"
-    "⚠️ Input neighbors are expected in CSR format. See demo.py for \n"
-    "examples of how to build such structure from typical scikit-learn\n"
-    " K-NN or radius-NN search outputs.\n\n"
-    "Parameters\n"
-    "----------\n"
-    ":param xyz: [n_points, 3] float32 2D array\n"
-    "    3D point coordinates\n"
-    ":param nn: [num_neighborhoods] uint32 1D array\n"
-    "    Flattened neighbor indices. Make sure those are all positive,\n"
-    "    '-1' indices will either crash or silently compute incorrect\n"
-    "    features\n"
-    ":param nn_ptr: [n_points+1] uint32 1D array\n"
-    "    Pointers wrt `nn`. More specifically, the neighbors of point `i`\n"
-    "    are `nn[nn_ptr[i]:nn_ptr[i + 1]]`\n"
-    ":param k_min: (optional, default=1) int\n"
-    "    Minimum number of neighbors to consider for features\n"
-    "    computation. If a point has less, it will be given 0 features\n"
-    ":param k_step: (optional, default=-1) int\n"
-    "    Step size to take when searching for the optimal neighborhood\n"
-    "    size for each point, following:\n"
-    "    http://lareg.ensg.eu/labos/matis/pdf/articles_revues/2015/isprs_wjhm_15.pdf\n"
-    "    If k_step < 1, pgeof will not search for the optimal\n"
-    "    neighborhood and features will be computed based on the all\n"
-    "    available neighbors for each point\n"
-    ":param k_min_search: (optional, default=10) int\n"
-    "    Minimum neighborhood size at which to start when searching for\n"
-    "    the optimal neighborhood size for each point. It is advised to\n"
-    "    use a value of 10 or higher, for geometric features robustness\n"
-    ":param verbose: (optional, default=False) bool\n"
-    "    Whether computation progress should be printed out\n";
-
-static PyMethodDef pgeof_methods[] = {
-    {"pgeof", (PyCFunction) pgeof_cpy, METH_VARARGS | METH_KEYWORDS, pgeof_doc},
-    {NULL, NULL, 0, NULL}
-};
-
-/* module initialization */
-
-static struct PyModuleDef pgeof_module = {
-    PyModuleDef_HEAD_INIT,
-    "pgeof", /* name of module */
-    "Pointwise geometric feature from point cloud", /* module documentation, may be null */
-    -1,   /* size of per-interpreter state of the module,
-             or -1 if the module keeps state in global variables. */
-    pgeof_methods, /* actual methods in the module */
-    NULL, /* multi-phase initialization, may be null */
-    NULL, /* traversal function, may be null */
-    NULL, /* clearing function, may be null */
-    NULL  /* freeing function, may be null */
-};
-
-PyMODINIT_FUNC
-PyInit_pgeof(void)
-{
-    import_array() /* IMPORTANT: this must be called to use numpy array */
-
-    PyObject* m;
-
-    /* create the module */
-    m = PyModule_Create(&pgeof_module);
-    if (!m){ return NULL; }
-
-    return m;
-}
diff --git a/src/pgeof_ext.cpp b/src/pgeof_ext.cpp
new file mode 100644
index 0000000..679c8f1
--- /dev/null
+++ b/src/pgeof_ext.cpp
@@ -0,0 +1,179 @@
+
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/pair.h>
+#include <nanobind/stl/vector.h>
+
+#include "nn_search.hpp"
+#include "pgeof.hpp"
+
+namespace nb = nanobind;
+using namespace nb::literals;
+
+NB_MODULE(pgeof_ext, m)
+{
+    m.doc() =
+        "Compute, for each point in a 3D point cloud, local geometric "
+        "features "
+        "in parallel on CPU";
+    nb::enum_<pgeof::EFeatureID>(m, "EFeatureID")
+        .value("Linearity", pgeof::EFeatureID::Linearity)
+        .value("Planarity", pgeof::EFeatureID::Planarity)
+        .value("Scattering", pgeof::EFeatureID::Scattering)
+        .value("VerticalityPGEOF", pgeof::EFeatureID::VerticalityPGEOF)
+        .value("Normal_x", pgeof::EFeatureID::Normal_x)
+        .value("Normal_y", pgeof::EFeatureID::Normal_y)
+        .value("Normal_z", pgeof::EFeatureID::Normal_z)
+        .value("Length", pgeof::EFeatureID::Length)
+        .value("Surface", pgeof::EFeatureID::Surface)
+        .value("Volume", pgeof::EFeatureID::Volume)
+        .value("Curvature", pgeof::EFeatureID::Curvature)
+        .value("K_optimal", pgeof::EFeatureID::K_optimal)  // TODO: remove or handle
+        .value("Verticality", pgeof::EFeatureID::Verticality)
+        .value("Eigentropy", pgeof::EFeatureID::Eigentropy)
+        .export_values();
+    m.def(
+        "compute_features", &pgeof::compute_geometric_features<float>, "xyz"_a.noconvert(), "nn"_a.noconvert(),
+        "nn_ptr"_a.noconvert(), "k_min"_a = 1, "verbose"_a = false, R"(
+            Compute a set of geometric features for a point cloud from a precomputed list of neighbors.
+
+            * The following features are computed:
+            - linearity
+            - planarity
+            - scattering
+            - verticality
+            - normal vector (oriented towards positive z-coordinates)
+            - length
+            - surface
+            - volume
+            - curvature
+            :param xyz: The point cloud. A numpy array of shape (n, 3).
+            :param nn: Integer 1D array. Flattened neighbor indices. Make sure those are all positive,
+            '-1' indices will either crash or silently compute incorrect features.
+            :param nn_ptr: [n_points+1] Integer 1D array. Pointers wrt 'nn'. More specifically, the neighbors of point 'i'
+            are 'nn[nn_ptr[i]:nn_ptr[i + 1]]'.
+            :param k_min: Minimum number of neighbors to consider for features computation. If a point has less,
+            its features will be a set of '0' values.
+            :param verbose: Whether computation progress should be printed out
+            :return: the geometric features associated with each point's neighborhood in a (num_points, features_count) numpy array.
+        )");
+    m.def(
+        "compute_features_multiscale", &pgeof::compute_geometric_features_multiscale<float>, "xyz"_a.noconvert(),
+        "nn"_a.noconvert(), "nn_ptr"_a.noconvert(), "k_scales"_a, "verbose"_a = false, R"(
+            Compute a set of geometric features for a point cloud in a multiscale fashion.
+            
+            * The following features are computed:
+            - linearity
+            - planarity
+            - scattering
+            - verticality
+            - normal vector (oriented towards positive z-coordinates)
+            - length
+            - surface
+            - volume
+            - curvature
+            
+            :param xyz: The point cloud. A numpy array of shape (n, 3).
+            :param nn: Integer 1D array. Flattened neighbor indices. Make sure those are all positive,
+            '-1' indices will either crash or silently compute incorrect features.
+            :param nn_ptr: [n_points+1] Integer 1D array. Pointers wrt 'nn'. More specifically, the neighbors of point 'i'
+            are 'nn[nn_ptr[i]:nn_ptr[i + 1]]'.
+            :param k_scale: Array of number of neighbors to consider for features computation. If a at a given scale, a point has
+            less features will be a set of '0' values.
+            :param verbose: Whether computation progress should be printed out
+            :return: Geometric features associated with each point's neighborhood in a (num_points, features_count, n_scales)
+            numpy array.
+        )");
+    m.def(
+        "compute_features_optimal", &pgeof::compute_geometric_features_optimal<float>, "xyz"_a.noconvert(),
+        "nn"_a.noconvert(), "nn_ptr"_a.noconvert(), "k_min"_a = 1, "k_step"_a = 1, "k_min_search"_a = 1,
+        "verbose"_a = false, R"(
+            Compute a set of geometric features for a point cloud using the optimal neighborhood selection described in
+            http://lareg.ensg.eu/labos/matis/pdf/articles_revues/2015/isprs_wjhm_15.pdf
+
+            * The following features are computed:
+            - linearity
+            - planarity
+            - scattering
+            - verticality
+            - normal vector (oriented towards positive z-coordinates)
+            - length
+            - surface
+            - volume
+            - curvature
+            - optimal_nn
+            :param xyz: the point cloud
+            :param nn: Integer 1D array. Flattened neighbor indices. Make sure those are all positive,
+            '-1' indices will either crash or silently compute incorrect features.
+            :param nn_ptr: [n_points+1] Integer 1D array. Pointers wrt 'nn'. More specifically, the neighbors of point 'i'
+            are 'nn[nn_ptr[i]:nn_ptr[i + 1]]'.
+            :param k_min: Minimum number of neighbors to consider for features computation. If a point has less,
+            its features will be a set of '0' values.
+            :param k_step: Step size to take when searching for the optimal neighborhood, size for each point following
+            Weinmann, 2015
+            :param k_min_search: Minimum neighborhood size at which to start when searching for the optimal neighborhood size for
+            each point. It is advised to use a value of 10 or higher, for geometric features robustness.
+            :param verbose: Whether computation progress should be printed out
+            :return: Geometric features associated with each point's neighborhood in a (num_points, features_count) numpy array.
+        )");
+    m.def("knn_search", &pgeof::nanoflann_knn_search<float>, "data"_a.noconvert(), "query"_a.noconvert(), "knn"_a, R"(
+        Given two point clouds, compute for each point present in one of the point cloud 
+        the N closest points in the other point cloud
+
+        It should be faster than scipy.spatial.KDTree for this task.
+        
+        :param data: the reference point cloud. A numpy array of shape (n, 3).
+        :param query: the point cloud used for the queries. A numpy array of shape (n, 3).
+        :param knn: the number of neighbors to take into account for each point.
+        :return: a pair of arrays, both of size (n_points x knn), the first one contains the indices of each neighbor, the
+        second one the square distances between the query point and each of its neighbors.
+    )");
+    m.def(
+        "radius_search", &pgeof::nanoflann_radius_search<float>, "data"_a.noconvert(), "query"_a.noconvert(),
+        "search_radius"_a, "knn"_a = 200, R"(
+            Search for the points within a specified sphere in a point cloud.
+            
+            It could be a fallback replacement for FRNN into SuperPointTransformer code base.
+            It should be faster than scipy.spatial.KDTree for this task.
+            
+            :param data: the reference point cloud. A numpy array of shape (n, 3).
+            :param query: the point cloud used for the queries (sphere centers). A numpy array of shape (n, 3).
+            :param search_radius: the search radius.
+            :param max_knn: the maximum number of neighbors to fetch inside the radius. The central point is included. Fixing a
+            reasonable max number of neighbors prevents running OOM for large radius/dense point clouds.
+            :return: a pair of arrays, both of size (n_points x knn), the first one contains the 'indices' of each neighbor,
+            the second one the 'square_distances' between the query point and each neighbor. Point having a number of neighbors <
+            'max_knn' inside the 'search_radius' will have their 'indices' and and 'square_distances' filled respectively with
+            '-1' and 'O' for any missing neighbor.
+        )");
+    m.def(
+        "compute_features_selected", &pgeof::compute_geometric_features_selected<double>, "xyz"_a.noconvert(),
+        "search_radius"_a, "knn"_a, "selected_features"_a, R"(
+            Compute a selected set of geometric features for a point cloud via radius search.
+
+            This function aims to mimick the behavior of jakteristics and provide an efficient way
+            to compute a limited set of features (double precision version).
+            
+            :param xyz: the point cloud. A numpy array of shape (n, 3).
+            :param search_radius: the search radius. A numpy array of shape (n, 3).
+            :param max_knn: the maximum number of neighbors to fetch inside the sphere. The central point is included. Fixing a
+            reasonable max number of neighbors prevents running OOM for large radius/dense point clouds.
+            :param selected_features: List of selected features. See EFeatureID
+            :return: Geometric features associated with each point's neighborhood in a (num_points, features_count) numpy array.
+        )");
+    m.def(
+        "compute_features_selected", &pgeof::compute_geometric_features_selected<float>, "xyz"_a.noconvert(),
+        "search_radius"_a, "knn"_a, "selected_features"_a, R"(
+            Compute a selected set of geometric features for a point cloud via radius search.
+
+            This function aims to mimic the behavior of jakteristics and provide an efficient way
+            to compute a limited set of features (float precision version).
+            
+            :param xyz: the point cloud
+            :param search_radius: the search radius.
+            :param max_knn: the maximum number of neighbors to fetch inside the sphere. The central point is included. Fixing a
+            reasonable max number of neighbors prevents running OOM for large radius/dense point clouds.
+            :param selected_features: List of selected features. See EFeatureID
+            :return: Geometric features associated with each point's neighborhood in a (num_points, features_count) numpy array.
+        )");
+}
\ No newline at end of file
diff --git a/third_party/eigen b/third_party/eigen
new file mode 160000
index 0000000..287c801
--- /dev/null
+++ b/third_party/eigen
@@ -0,0 +1,1 @@
+Subproject commit 287c8017808fb4c7a651ba4c02363a773e2f0c46
diff --git a/third_party/nanoflann b/third_party/nanoflann
new file mode 160000
index 0000000..923c2ac
--- /dev/null
+++ b/third_party/nanoflann
@@ -0,0 +1,1 @@
+Subproject commit 923c2ac16a955317054c296d298425d33c0cb9de
diff --git a/third_party/taskflow b/third_party/taskflow
new file mode 160000
index 0000000..12f8bd4
--- /dev/null
+++ b/third_party/taskflow
@@ -0,0 +1,1 @@
+Subproject commit 12f8bd4e970ab27fd3dee3bffa24b5b48b54ba39


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 0000000..64f466a
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,29 @@
+name: Tests
+
+on:
+  - push
+  - pull_request
+
+jobs:
+  build:
+    name: Tests on ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-22.04, windows-2022] 
+        python-version: ['3.9', '3.10', '3.11', '3.12']
+    steps:
+    - uses: actions/checkout@v4
+      with:
+        submodules: True
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v4
+      with:
+        python-version: ${{ matrix.python-version }}
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        python -m pip install tox tox-gh-actions
+    - name: Test with tox
+      run: tox
\ No newline at end of file
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/bench_jakteristics.py b/tests/bench_jakteristics.py
new file mode 100644
index 0000000..298b1c8
--- /dev/null
+++ b/tests/bench_jakteristics.py
@@ -0,0 +1,45 @@
+import numpy as np
+import pytest
+
+import pgeof
+from pgeof import EFeatureID
+
+# Skip if jakteristics import fail
+# it should fail on darwin (macOS) systems
+jakteristics = pytest.importorskip("jakteristics")
+
+
+@pytest.fixture
+def random_point_cloud():
+    return np.random.rand(1000000, 3) * 100
+
+
+@pytest.mark.benchmark(group="feature-computation-jak", disable_gc=True, warmup=True)
+def test_bench_jak(benchmark, random_point_cloud):
+    knn = 50
+    dist = 5.0
+
+    def _to_bench_feat():
+        _ = jakteristics.compute_features(
+            random_point_cloud,
+            dist,
+            kdtree=None,
+            num_threads=-1,
+            max_k_neighbors=knn,
+            feature_names=["verticality"],
+        )
+
+    benchmark(_to_bench_feat)
+
+
+@pytest.mark.benchmark(group="feature-computation-jak", disable_gc=True, warmup=True)
+def test_pgeof(benchmark, random_point_cloud):
+    knn = 50
+    dist = 5.0
+
+    def _to_bench_feat():
+        _ = pgeof.compute_features_selected(
+            random_point_cloud, dist, knn, [EFeatureID.Verticality]
+        )
+
+    benchmark(_to_bench_feat)
diff --git a/tests/bench_knn.py b/tests/bench_knn.py
new file mode 100644
index 0000000..6aa4c83
--- /dev/null
+++ b/tests/bench_knn.py
@@ -0,0 +1,54 @@
+import numpy as np
+import pytest
+from scipy.spatial import KDTree
+
+import pgeof
+
+
+@pytest.fixture
+def random_point_cloud():
+    return np.random.rand(10000000, 3).astype("float32")
+
+
+@pytest.mark.benchmark(group="knn", disable_gc=True, warmup=True)
+def test_knn_scipy(benchmark, random_point_cloud):
+    knn = 50
+
+    def _to_bench():
+        tree = KDTree(random_point_cloud)
+        _ = tree.query(random_point_cloud, k=knn, workers=-1)
+
+    benchmark(_to_bench)
+
+
+@pytest.mark.benchmark(group="knn", disable_gc=True, warmup=True)
+def test_knn_pgeof(benchmark, random_point_cloud):
+    knn = 50
+
+    def _to_bench():
+        _ = pgeof.knn_search(random_point_cloud, random_point_cloud, knn)
+
+    benchmark(_to_bench)
+
+
+@pytest.mark.benchmark(group="radius-search", disable_gc=True, warmup=True)
+def test_radius_scipy(benchmark, random_point_cloud):
+    max_knn = 30
+    radius = 0.2
+
+    def _to_bench():
+        tree = KDTree(random_point_cloud)
+        _ = tree.query(random_point_cloud, k=max_knn, distance_upper_bound=radius, workers=-1)
+
+    benchmark(_to_bench)
+
+
+@pytest.mark.benchmark(group="radius-search", disable_gc=True, warmup=True)
+def test_radius_pgeof(benchmark, random_point_cloud):
+    max_knn = 30
+    radius = 0.2
+
+    def _to_bench():
+        _ = pgeof.radius_search(random_point_cloud, random_point_cloud, radius, max_knn)
+
+    benchmark(_to_bench)
diff --git a/tests/helpers.py b/tests/helpers.py
new file mode 100644
index 0000000..fce33cb
--- /dev/null
+++ b/tests/helpers.py
@@ -0,0 +1,23 @@
+import numpy as np
+from scipy.spatial import KDTree
+
+
+def random_nn(num_points, k):
+    # Generate a random synthetic point cloud
+    xyz = np.random.rand(num_points, 3)
+
+    # Converting k-nearest neighbors to CSR format
+    kneigh = KDTree(xyz).query(xyz, k=k, workers=-1)
+    nn_ptr = np.arange(num_points + 1) * k
+    nn = kneigh[1].flatten()
+
+    # Make sure xyz are float32 and nn and nn_ptr are uint32
+    xyz = xyz.astype("float32")
+    nn_ptr = nn_ptr.astype("uint32")
+    nn = nn.astype("uint32")
+
+    # Make sure arrays are contiguous (C-order) and not Fortran-order
+    xyz = np.ascontiguousarray(xyz)
+    nn_ptr = np.ascontiguousarray(nn_ptr)
+    nn = np.ascontiguousarray(nn)
+    return xyz, nn, nn_ptr
\ No newline at end of file
diff --git a/tests/test_pgeof.py b/tests/test_pgeof.py
new file mode 100644
index 0000000..7610df4
--- /dev/null
+++ b/tests/test_pgeof.py
@@ -0,0 +1,44 @@
+import numpy as np
+from scipy.spatial import KDTree
+
+import pgeof
+from tests.helpers import random_nn
+
+
+def test_knn():
+    knn = 10
+    xyz = np.random.rand(1000, 3)
+    xyz = xyz.astype("float32")
+    tree = KDTree(xyz)
+    _, k_legacy = tree.query(xyz, k=knn, workers=-1)
+    k_new, _ = pgeof.knn_search(xyz, xyz, knn)
+    np.testing.assert_equal(k_legacy, k_new)
+
+def test_radius_search():
+    knn = 10
+    radius = 0.2
+    xyz = np.random.rand(1000, 3)
+    xyz = xyz.astype("float32")
+    tree = KDTree(xyz)
+    _, k_legacy = tree.query(xyz, k=knn, distance_upper_bound=radius, workers=-1)
+    k_legacy[k_legacy == xyz.shape[0]] = -1
+    k_new, _ = pgeof.radius_search(xyz, xyz, radius, knn)
+    np.testing.assert_equal(k_legacy, k_new)
+
+def test_pgeof_multiscale():
+    # Generate a random synthetic point cloud and NNs
+    xyz, nn, nn_ptr = random_nn(10000, 50)
+
+    # with pytest.raises(ValueError):
+    # scales = np.array(
+    #   [20, 50]
+    # )  # scales in decreasing order in order to raise the exception
+    # multi = pgeof.compute_features_multiscale(xyz, nn, nn_ptr, scales, False)
+    scales = np.array(
+        [50, 20]
+    )  # scales in decreasing order in order to raise the exception
+    multi = pgeof.compute_features_multiscale(xyz, nn, nn_ptr, np.flip(scales), False)
+    simple = pgeof.compute_features(xyz, nn, nn_ptr, 50, False)
+    multi_simple = pgeof.compute_features_multiscale(xyz, nn, nn_ptr, [20], False)
+    np.testing.assert_allclose(multi[:, 0], multi_simple[:, 0], 1e-1, 1e-5)
+    np.testing.assert_allclose(multi[:, 1], simple, 1e-1, 1e-5)
\ No newline at end of file