Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "Blobs.h"
       4             : #include "BSplines.h"
       5             : #include "Logger.h"
       6             : #include "Timer.h"
       7             : 
       8             : #include <array>
       9             : 
      10             : namespace elsa
      11             : {
      12             :     namespace detail
      13             :     {
      14             :         template <typename data_t, index_t N>
      15           0 :         constexpr std::array<data_t, N> blob_lut(ProjectedBlob<data_t> blob)
      16             :         {
      17           0 :             Logger::get("blob_lut")->debug("Calculating lut");
      18             : 
      19             :             std::array<data_t, N> lut;
      20             : 
      21           0 :             auto t = static_cast<data_t>(0);
      22           0 :             const auto step = blob.radius() / N;
      23             : 
      24           0 :             for (std::size_t i = 0; i < N; ++i) {
      25           0 :                 lut[i] = blob(t);
      26           0 :                 t += step;
      27             :             }
      28             : 
      29           0 :             return lut;
      30             :         }
      31             : 
      32             :         template <typename data_t, index_t N>
      33           0 :         constexpr std::array<data_t, N> bspline_lut(ProjectedBSpline<data_t> bspline)
      34             :         {
      35           0 :             Logger::get("bspline_lut")->debug("Calculating lut");
      36             : 
      37             :             std::array<data_t, N> lut;
      38             : 
      39           0 :             auto t = static_cast<data_t>(0);
      40           0 :             const auto step = 2. / N;
      41             : 
      42           0 :             for (std::size_t i = 0; i < N; ++i) {
      43           0 :                 lut[i] = bspline(t);
      44           0 :                 t += step;
      45             :             }
      46             : 
      47           0 :             return lut;
      48             :         }
      49             : 
      50             :         template <typename data_t>
      51           0 :         data_t lerp(data_t a, SelfType_t<data_t> b, SelfType_t<data_t> t)
      52             :         {
      53           0 :             if ((a <= 0 && b >= 0) || (a >= 0 && b <= 0))
      54           0 :                 return t * b + (1 - t) * a;
      55             : 
      56           0 :             if (t == 1)
      57           0 :                 return b;
      58             : 
      59           0 :             const data_t x = a + t * (b - a);
      60             : 
      61           0 :             if ((t > 1) == (b > a))
      62           0 :                 return b < x ? x : b;
      63             :             else
      64           0 :                 return x < b ? x : b;
      65             :         }
      66             :     } // namespace detail
      67             : 
      68             :     template <typename data_t, std::size_t N>
      69             :     class Lut
      70             :     {
      71             :     public:
      72           0 :         constexpr Lut(std::array<data_t, N> data) : data_(std::move(data)) {}
      73             : 
      74             :         template <typename T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
      75           0 :         constexpr data_t operator()(T index) const
      76             :         {
      77           0 :             if (index < 0 || index > N) {
      78           0 :                 return 0;
      79             :             }
      80             : 
      81           0 :             return data_[index];
      82             :         }
      83             : 
      84             :         /// TODO: Handle boundary conditions
      85             :         /// lerp(last, last+1, t), for some t > 0, yields f(last) / 2, as f(last + 1) = 0,
      86             :         /// this should be handled
      87             :         template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
      88           0 :         constexpr data_t operator()(T index) const
      89             :         {
      90           0 :             if (index < 0 || index > N) {
      91           0 :                 return 0;
      92             :             }
      93             : 
      94             :             // Get the two closes indices
      95           0 :             const auto a = static_cast<std::size_t>(std::floor(index));
      96           0 :             const auto b = static_cast<std::size_t>(std::ceil(index));
      97             : 
      98             :             // Get the function values
      99           0 :             const auto fa = a == N ? 0 : data_[a];
     100           0 :             const auto fb = b == N ? 0 : data_[b];
     101             : 
     102           0 :             auto ret = detail::lerp(fa, fb, index - static_cast<data_t>(a));
     103             : 
     104             :             // Bilinear interpolation
     105           0 :             return detail::lerp(fa, fb, index - static_cast<data_t>(a));
     106             :         }
     107             : 
     108             :     private:
     109             :         std::array<data_t, N> data_;
     110             :     };
     111             : 
     112             :     // User defined deduction guide
     113             :     template <typename data_t, std::size_t N>
     114             :     Lut(std::array<data_t, N>) -> Lut<data_t, N>;
     115             : 
     116             :     template <typename data_t, index_t N>
     117             :     class ProjectedBlobLut
     118             :     {
     119             :     public:
     120           0 :         constexpr ProjectedBlobLut(data_t radius, SelfType_t<data_t> alpha,
     121             :                                    SelfType_t<data_t> order)
     122           0 :             : blob_(radius, alpha, order), lut_(detail::blob_lut<data_t, N>(blob_))
     123             :         {
     124           0 :         }
     125             : 
     126           0 :         constexpr data_t radius() const { return blob_.radius(); }
     127             : 
     128             :         constexpr data_t alpha() const { return blob_.alpha(); }
     129             : 
     130             :         constexpr data_t order() const { return blob_.order(); }
     131             : 
     132           0 :         constexpr data_t operator()(data_t distance) const
     133             :         {
     134           0 :             return lut_((distance / blob_.radius()) * N);
     135             :         }
     136             : 
     137             :     private:
     138             :         ProjectedBlob<data_t> blob_;
     139             :         Lut<data_t, N> lut_;
     140             :     };
     141             : 
     142             :     template <typename data_t, index_t N>
     143             :     class ProjectedBSplineLut
     144             :     {
     145             :     public:
     146           0 :         constexpr ProjectedBSplineLut(int dim, int degree)
     147           0 :             : bspline_(dim, degree), lut_(detail::bspline_lut<data_t, N>(bspline_))
     148             :         {
     149           0 :         }
     150             : 
     151             :         constexpr data_t order() const { return bspline_.order(); }
     152             : 
     153           0 :         constexpr data_t operator()(data_t distance) const
     154             :         {
     155           0 :             return lut_((std::abs(distance) / 2.) * N);
     156             :         }
     157             : 
     158             :     private:
     159             :         ProjectedBSpline<data_t> bspline_;
     160             :         Lut<data_t, N> lut_;
     161             :     };
     162             : } // namespace elsa