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             :         constexpr std::array<data_t, N> blob_lut(ProjectedBlob<data_t> blob)
      16         875 :         {
      17         875 :             Logger::get("blob_lut")->debug("Calculating lut");
      18             : 
      19         875 :             std::array<data_t, N> lut;
      20             : 
      21         875 :             auto t = static_cast<data_t>(0);
      22         875 :             const auto step = blob.radius() / N;
      23             : 
      24       88375 :             for (std::size_t i = 0; i < N; ++i) {
      25       87500 :                 lut[i] = blob(t);
      26       87500 :                 t += step;
      27       87500 :             }
      28             : 
      29         875 :             return lut;
      30         875 :         }
      31             : 
      32             :         template <typename data_t, index_t N>
      33             :         constexpr std::array<data_t, N> bspline_lut(ProjectedBSpline<data_t> bspline)
      34          16 :         {
      35          16 :             Logger::get("bspline_lut")->debug("Calculating lut");
      36             : 
      37          16 :             std::array<data_t, N> lut;
      38             : 
      39          16 :             auto t = static_cast<data_t>(0);
      40          16 :             const auto step = 2. / N;
      41             : 
      42         816 :             for (std::size_t i = 0; i < N; ++i) {
      43         800 :                 lut[i] = bspline(t);
      44         800 :                 t += step;
      45         800 :             }
      46             : 
      47          16 :             return lut;
      48          16 :         }
      49             : 
      50             :         template <typename data_t>
      51             :         data_t lerp(data_t a, SelfType_t<data_t> b, SelfType_t<data_t> t)
      52       84344 :         {
      53       84344 :             if ((a <= 0 && b >= 0) || (a >= 0 && b <= 0))
      54        3078 :                 return t * b + (1 - t) * a;
      55             : 
      56       81266 :             if (t == 1)
      57          32 :                 return b;
      58             : 
      59       81234 :             const data_t x = a + t * (b - a);
      60             : 
      61       81234 :             if ((t > 1) == (b > a))
      62       79611 :                 return b < x ? x : b;
      63        1623 :             else
      64        1623 :                 return x < b ? x : b;
      65       81234 :         }
      66             :     } // namespace detail
      67             : 
      68             :     template <typename data_t, std::size_t N>
      69             :     class Lut
      70             :     {
      71             :     public:
      72         905 :         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             :         constexpr data_t operator()(T index) const
      76         400 :         {
      77         400 :             if (index < 0 || index > N) {
      78           0 :                 return 0;
      79           0 :             }
      80             : 
      81         400 :             return data_[index];
      82         400 :         }
      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             :         constexpr data_t operator()(T index) const
      89       55728 :         {
      90       55740 :             if (index < 0 || index > N) {
      91       13432 :                 return 0;
      92       13432 :             }
      93             : 
      94             :             // Get the two closes indices
      95       42296 :             const auto a = static_cast<std::size_t>(std::floor(index));
      96       42296 :             const auto b = static_cast<std::size_t>(std::ceil(index));
      97             : 
      98             :             // Get the function values
      99       42296 :             const auto fa = a == N ? 0 : data_[a];
     100       42296 :             const auto fb = b == N ? 0 : data_[b];
     101             : 
     102       42296 :             auto ret = detail::lerp(fa, fb, index - static_cast<data_t>(a));
     103             : 
     104             :             // Bilinear interpolation
     105       42296 :             return detail::lerp(fa, fb, index - static_cast<data_t>(a));
     106       42296 :         }
     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             :         constexpr ProjectedBlobLut(data_t radius, SelfType_t<data_t> alpha,
     121             :                                    SelfType_t<data_t> order)
     122             :             : blob_(radius, alpha, order), lut_(detail::blob_lut<data_t, N>(blob_))
     123         875 :         {
     124         875 :         }
     125             : 
     126        1611 :         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             :         constexpr data_t operator()(data_t distance) const
     133       51603 :         {
     134       51603 :             return lut_((distance / blob_.radius()) * N);
     135       51603 :         }
     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             :         constexpr ProjectedBSplineLut(int dim, int degree)
     147             :             : bspline_(dim, degree), lut_(detail::bspline_lut<data_t, N>(bspline_))
     148          16 :         {
     149          16 :         }
     150             : 
     151             :         constexpr data_t order() const { return bspline_.order(); }
     152             : 
     153             :         constexpr data_t operator()(data_t distance) const
     154        3200 :         {
     155        3200 :             return lut_((std::abs(distance) / 2.) * N);
     156        3200 :         }
     157             : 
     158             :     private:
     159             :         ProjectedBSpline<data_t> bspline_;
     160             :         Lut<data_t, N> lut_;
     161             :     };
     162             : } // namespace elsa