metrics.hpp

Go to the documentation of this file.

/* Copyright (c) 2025 Otto Link. Distributed under the terms of the GNU General
   Public License. The full license is in the file LICENSE, distributed with
   this software. */
#pragma once
#include <cstddef>
#include <optional>
 
#include "point_sampler/internal/nanoflann_adaptator.hpp"
#include "point_sampler/point.hpp"
 
namespace ps
{
 
template <typename T, size_t N>
std::pair<std::vector<T>, std::vector<T>> angle_distribution_neighbors(
    const std::vector<Point<T, N>> &points,
    T                               bin_width,
    size_t                          k_neighbors = 8)
{
  if constexpr (N < 2)
    throw std::runtime_error("Angle distribution requires dimension >= 2");
 
  size_t         num_bins = static_cast<size_t>(std::ceil(M_PI / bin_width));
  std::vector<T> angles(num_bins);
  std::vector<T> g_theta(num_bins, T(0));
 
  // Precompute bin centers
  for (size_t i = 0; i < num_bins; ++i)
    angles[i] = (i + T(0.5)) * bin_width;
 
  // Get nearest neighbors once
  auto neighbors = nearest_neighbors_indices(points, k_neighbors);
 
  // Loop over all points
  for (size_t i = 0; i < points.size(); ++i)
  {
    const auto &nbrs = neighbors[i];
    size_t      n_nbrs = nbrs.size();
 
    // Compute angles between all neighbor pairs
    for (size_t a = 0; a < n_nbrs; ++a)
    {
      for (size_t b = a + 1; b < n_nbrs; ++b)
      {
        size_t j = nbrs[a];
        size_t k = nbrs[b];
 
        // Vectors from center i to j and k
        std::array<T, N> v1, v2;
        for (size_t d = 0; d < N; ++d)
        {
          v1[d] = points[j][d] - points[i][d];
          v2[d] = points[k][d] - points[i][d];
        }
 
        // Compute angle
        T dot = T(0), norm1 = T(0), norm2 = T(0);
        for (size_t d = 0; d < N; ++d)
        {
          dot += v1[d] * v2[d];
          norm1 += v1[d] * v1[d];
          norm2 += v2[d] * v2[d];
        }
 
        if (norm1 > T(0) && norm2 > T(0))
        {
          T cos_theta = std::clamp(dot / (std::sqrt(norm1) * std::sqrt(norm2)),
                                   T(-1),
                                   T(1));
          T theta = std::acos(cos_theta);
 
          size_t bin = static_cast<size_t>(theta / bin_width);
          if (bin < num_bins)
            g_theta[bin] += T(1);
        }
      }
    }
  }
 
  // Normalize histogram
  T total = std::accumulate(g_theta.begin(), g_theta.end(), T(0));
  if (total > T(0))
    for (auto &val : g_theta)
      val /= total;
 
  return {angles, g_theta};
}
 
template <typename T, size_t N>
std::vector<T> distance_to_boundary(const std::vector<Point<T, N>>       &points,
                                    const std::array<std::pair<T, T>, N> &axis_ranges)
{
  std::vector<T> distances;
  distances.reserve(points.size());
 
  for (const auto &p : points)
  {
    // Find min distance to a boundary plane
    T min_dist = std::numeric_limits<T>::max();
 
    for (size_t d = 0; d < N; ++d)
    {
      T dist_to_min = std::abs(p[d] - axis_ranges[d].first);
      T dist_to_max = std::abs(axis_ranges[d].second - p[d]);
 
      min_dist = std::min({min_dist, dist_to_min, dist_to_max});
    }
 
    distances.push_back(min_dist);
  }
 
  return distances;
}
 
template <typename T, size_t N>
std::vector<T> first_neighbor_distance_squared(std::vector<Point<T, N>> &points)
{
  PointCloudAdaptor<T, N> adaptor(points);
  KDTree<T, N>            index(N, adaptor);
  index.buildIndex();
 
  std::vector<T> distance_sq;
  distance_sq.reserve(points.size());
 
  // first neighbor search only
  const size_t k_neighbors = 1;
 
  for (size_t i = 0; i < points.size(); ++i)
  {
    const auto &p = points[i];
 
    std::vector<size_t> ret_indexes(k_neighbors + 1);
    std::vector<T>      out_dists_sqr(k_neighbors + 1);
 
    nanoflann::KNNResultSet<T> result_set(k_neighbors + 1);
    result_set.init(ret_indexes.data(), out_dists_sqr.data());
 
    std::array<T, N> query;
    for (size_t d = 0; d < N; ++d)
      query[d] = p[d];
 
    index.findNeighbors(result_set, query.data(), nanoflann::SearchParameters());
 
    T dist_sq = 0.f;
    if (ret_indexes.size() > 1)
    {
      const auto &q = points[ret_indexes[1]];
      auto        delta = p - q;
      dist_sq = length_squared(delta) + T(1e-6);
    }
 
    distance_sq.push_back(dist_sq);
  }
 
  return distance_sq;
}
 
template <typename T, size_t N>
std::vector<T> local_density_knn(const std::vector<Point<T, N>> &points, size_t k = 8)
{
  auto neighbors = nearest_neighbors_indices(points, k);
 
  std::vector<T> densities(points.size());
 
  // Compute unit N-ball volume using gamma function
  T volume_unit_ball = std::pow(M_PI, T(N) / 2) / std::tgamma(T(N) / 2 + 1);
 
  for (size_t i = 0; i < points.size(); ++i)
  {
    // Distance to k-th nearest neighbor
    const auto &nk_idx = neighbors[i].back();
    T           dist2 = 0;
    for (size_t d = 0; d < N; ++d)
    {
      T diff = points[i][d] - points[nk_idx][d];
      dist2 += diff * diff;
    }
    T r = std::sqrt(dist2);
 
    // Density = k / (V_N * r^N)
    T volume = volume_unit_ball * std::pow(r, T(N));
    densities[i] = k / volume;
  }
 
  return densities;
}
 
template <typename T, size_t N>
std::vector<std::vector<size_t>> nearest_neighbors_indices(
    const std::vector<Point<T, N>> &points,
    size_t                          k_neighbors = 8)
{
  PointCloudAdaptor<T, N> adaptor(points);
  KDTree<T, N>            index(N, adaptor);
  index.buildIndex();
 
  std::vector<std::vector<size_t>> all_neighbors;
  all_neighbors.resize(points.size());
 
  for (size_t i = 0; i < points.size(); ++i)
  {
    const auto &p = points[i];
 
    std::vector<size_t> ret_indexes(k_neighbors + 1);
    std::vector<T>      out_dists_sqr(k_neighbors + 1);
 
    nanoflann::KNNResultSet<T> result_set(k_neighbors + 1);
    result_set.init(ret_indexes.data(), out_dists_sqr.data());
 
    std::array<T, N> query;
    for (size_t d = 0; d < N; ++d)
      query[d] = p[d];
 
    index.findNeighbors(result_set, query.data(), nanoflann::SearchParameters());
 
    // Skip self at index 0
    all_neighbors[i].assign(ret_indexes.begin() + 1, ret_indexes.end());
  }
 
  return all_neighbors;
}
 
template <typename T, size_t N>
std::pair<std::vector<T>, std::vector<T>> radial_distribution(
    const std::vector<Point<T, N>>       &points,
    const std::array<std::pair<T, T>, N> &axis_ranges,
    T                                     bin_width,
    T                                     max_distance)
{
  size_t         num_bins = static_cast<size_t>(std::ceil(max_distance / bin_width));
  std::vector<T> radii(num_bins);
  std::vector<T> g(num_bins, T(0));
 
  // Precompute bin centers
  for (size_t i = 0; i < num_bins; ++i)
    radii[i] = (i + T(0.5)) * bin_width;
 
  // Compute volume of domain
  T volume = T(1);
  for (const auto &range : axis_ranges)
    volume *= (range.second - range.first);
 
  size_t n_points = points.size();
  T      density = static_cast<T>(n_points) / volume;
 
  // Histogram of pair distances
  for (size_t i = 0; i < n_points; ++i)
  {
    for (size_t j = i + 1; j < n_points; ++j)
    {
      T dist = distance(points[i], points[j]);
      if (dist < max_distance)
      {
        size_t bin = static_cast<size_t>(dist / bin_width);
        if (bin < num_bins)
          g[bin] += T(2); // count both i→j and
                          // j→i
      }
    }
  }
 
  // Normalize
  T norm_factor = (T)n_points * density;
 
  // Volume of spherical shell between r1 and r2 in N dimensions
  auto sphere_volume = [](T radius)
  { return std::pow(M_PI, N / 2.0) / std::tgamma(N / 2.0 + 1.0) * std::pow(radius, N); };
 
  for (size_t i = 0; i < num_bins; ++i)
  {
    T r_inner = i * bin_width;
    T r_outer = (i + 1) * bin_width;
 
    // Volume of spherical shell
    T shell_vol = sphere_volume(r_outer) - sphere_volume(r_inner);
    g[i] /= norm_factor * shell_vol;
  }
 
  return {radii, g};
}
 
} // namespace ps