Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "elsaDefines.h"
       4             : #include "Assertions.h"
       5             : #include "Math.hpp"
       6             : 
       7             : #include <cmath>
       8             : #include "spdlog/fmt/bundled/core.h"
       9             : #include "spdlog/fmt/fmt.h"
      10             : 
      11             : namespace elsa
      12             : {
      13             :     namespace bspline
      14             :     {
      15             :         /// @brief Evaluate the 1-dimensional B-Spline of degree m, with m + 2 equally spaced knots
      16             :         ///
      17             :         /// This is based on the implementation given as in e.g.:
      18             :         /// Fast B-spline transforms for continuous image representation and interpolation, Unser
      19             :         /// et. al. (1991), equation 2.2
      20             :         ///
      21             :         /// @param x coordinate to evaluate 1-dimensional B-Spline at
      22             :         /// @param m order of B-Spline
      23             :         template <typename data_t>
      24             :         constexpr data_t bspline1d_evaluate(data_t x, index_t m) noexcept
      25        6126 :         {
      26        6126 :             const auto xs = x + (m + 1) / data_t{2};
      27        6126 :             data_t res = 0;
      28             : 
      29       27630 :             for (int n = 0; n <= m + 1; ++n) {
      30       21504 :                 const auto tmp1 = math::heaviside<data_t>(xs - n, 0);
      31       21504 :                 const auto tmp2 = std::pow<data_t>(xs - n, m);
      32       21504 :                 const auto tmp3 = math::binom(m + 1, n);
      33       21504 :                 const auto tmp4 = std::pow<data_t>(-1, n) / math::factorial(m);
      34       21504 :                 res += tmp1 * tmp2 * tmp3 * tmp4;
      35       21504 :             }
      36             : 
      37        6126 :             return res;
      38        6126 :         }
      39             : 
      40             :         /// @brief Evaluate n-dimensional B-Spline of degree m. As B-Splines are separable, this is
      41             :         /// just the product of 1-dimensional B-Splines.
      42             :         ///
      43             :         /// @param x n-dimensional coordinate to evaluate B-Spline at
      44             :         /// @param m order of B-Spline
      45             :         template <typename data_t>
      46             :         constexpr data_t nd_bspline_evaluate(const Vector_t<data_t>& x, index_t m) noexcept
      47           0 :         {
      48           0 :             data_t res = bspline1d_evaluate(x[0], m);
      49           0 :             for (int i = 1; i < x.size(); ++i) {
      50           0 :                 res *= bspline1d_evaluate(x[i], m);
      51           0 :             }
      52           0 :             return res;
      53           0 :         }
      54             : 
      55             :         /// @brief Evaluate n-dimensional B-Spline at a given coordinate for the first dimension,
      56             :         /// and at the center (i.e. `0`) at all the other dimensions. This is particular useful as
      57             :         /// an approximation during the calculation of the line integral
      58             :         ///
      59             :         /// @param x 1-dimensional coordinate to evaluate B-Spline at
      60             :         /// @param m order of B-Spline
      61             :         /// @param dim dimension of B-Spline
      62             :         template <typename data_t>
      63             :         constexpr data_t nd_bspline_centered(data_t x, int m, int dim) noexcept
      64        4042 :         {
      65        4042 :             data_t res = bspline1d_evaluate<data_t>(x, m);
      66        6063 :             for (int i = 1; i < dim; ++i) {
      67        2021 :                 const auto inc = bspline1d_evaluate<data_t>(0., m);
      68        2021 :                 res *= inc;
      69        2021 :             }
      70        4042 :             return res;
      71        4042 :         }
      72             :     } // namespace bspline
      73             : 
      74             :     /// @brief Represent a B-Spline basis function of a given dimension and order
      75             :     template <typename data_t>
      76             :     class BSpline
      77             :     {
      78             :     public:
      79             :         BSpline(index_t dim, index_t order);
      80             : 
      81             :         data_t operator()(Vector_t<data_t> x);
      82             : 
      83             :         index_t order() const;
      84             : 
      85             :         index_t dim() const;
      86             : 
      87             :     private:
      88             :         /// Dimension of B-Spline
      89             :         index_t dim_;
      90             : 
      91             :         /// Order of B-Spline
      92             :         index_t order_;
      93             :     };
      94             : 
      95             :     /// @brief Represent a projected B-Spline basis function of a given dimension and order.
      96             :     /// Projected B-Splines are again B-Spline of n-1 dimensions. Using the fact, that B-Splines
      97             :     /// are close to symmetrical, we can approximate the projection only based on distance.
      98             :     template <typename data_t>
      99             :     class ProjectedBSpline
     100             :     {
     101             :     public:
     102             :         ProjectedBSpline(index_t dim, index_t order);
     103             : 
     104             :         data_t operator()(data_t s);
     105             : 
     106             :         index_t order() const;
     107             : 
     108             :         index_t dim() const;
     109             : 
     110             :     private:
     111             :         /// Dimension of B-Spline
     112             :         index_t dim_;
     113             : 
     114             :         /// Order of B-Spline
     115             :         index_t order_;
     116             :     };
     117             : } // namespace elsa