Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "Error.h"
       4             : #include "elsaDefines.h"
       5             : #include "Bessel.h"
       6             : 
       7             : namespace elsa
       8             : {
       9             :     template <class T>
      10             :     struct SelfType {
      11             :         using type = T;
      12             :     };
      13             : 
      14             :     template <class T>
      15             :     using SelfType_t = typename SelfType<T>::type;
      16             : 
      17             :     namespace blobs
      18             :     {
      19             :         // The devinition given by Lewitt (1992):
      20             :         // \f$
      21             :         // b_{m, alpha, a}(r) = \frac{I_m(\alpha * \sqrt{1 - (r/a)^2)}{I_m(\alpha)} * (\sqrt{1 -
      22             :         // (r/a)^2})^m \f$ for any \f$ 0 <= r <= a \f$.
      23             :         template <typename data_t>
      24             :         constexpr data_t blob_evaluate(data_t r, SelfType_t<data_t> a, SelfType_t<data_t> alpha,
      25             :                                        int m) noexcept
      26      368328 :         {
      27      368328 :             const auto w = static_cast<data_t>(1) - std::pow(r / a, static_cast<data_t>(2));
      28      368328 :             if (w >= 0) {
      29      367828 :                 const data_t Im1 = math::bessi(m, alpha);
      30      367828 :                 const data_t arg = std::sqrt(w);
      31      367828 :                 const data_t Im2 = math::bessi(m, alpha * arg);
      32             : 
      33      367828 :                 return Im2 / Im1 * std::pow(arg, m);
      34      367828 :             }
      35         500 :             return 0;
      36         500 :         }
      37             : 
      38             :         /// @brief Compute line integral of blob through a straight line
      39             :         /// The exact formulations are quite easily derived using the fact, that the line
      40             :         /// integral for a single blob is given by:
      41             :         ///
      42             :         /// \f$
      43             :         /// \frac{a}{I_m(\alpha)} \sqrt{\frac{2 \pi}{\alpha}} \left(\sqrt{1 -
      44             :         /// \left(\frac{s}{a}\right)^2}\right)^{m + 0.5} I_{m + 0.5}\left(\alpha \sqrt{1 -
      45             :         /// \left(\frac{s}{a}\right)^2}\right)
      46             :         /// \f$
      47             :         ///
      48             :         /// Then for a given order substitute I_{m + 0.5}(x) with the elementary function
      49             :         /// representations, as they can be found in Abramowitz and Stegun's Handbook of
      50             :         /// mathematical functions (1972):
      51             :         /// \f$ I_{0.5}(x) = \sqrt{\frac{2}{\pi x}} \sinh(x)\$f
      52             :         /// \f$ I_{1.5}(x) = \sqrt{\frac{2}{\pi x}} \left( \cosh(x) \frac{\sinh(x)}{x} \right) \$f
      53             :         /// \f$ I_{2.5}(x) = \sqrt{\frac{2}{\pi x}} \left(\left(\frac{3}{x^2} +
      54             :         /// \frac{1}{x}\right)\sinh(x) - \frac{3}{x^2} \cosh(x)\right)\$f
      55             :         ///
      56             :         /// Which will result in the below formulations. In theory using the recurrent relations,
      57             :         /// this could be extended to further orders, but is deemed unnecessary for now.
      58             :         ///
      59             :         /// TODO: What about alpha = 0? Check if you find anything in the literature for that.
      60             :         ///
      61             :         /// @param distance distance of blob center to straight line, in literature often referred
      62             :         /// to as `r`
      63             :         /// @param radius radius of blob, often referred to as `a`
      64             :         /// @param alpha smoothness factor of blob, expected to be larger than 0
      65             :         /// @param order order of Bessel function, often referred to as `m`
      66             :         ///
      67             :         /// Ref:
      68             :         /// Distance-Driven Projection and Backprojection for Spherically Symmetric Basis Functions
      69             :         /// in CT - Levakhina
      70             :         /// Spherically symmetric volume elements as basis functions for image reconstructions in
      71             :         /// computed laminography - P. Trampert
      72             :         /// Semi-Discrete Iteration Methods in X-Ray Tomography - Jonas Vogelgesang
      73             :         template <typename data_t>
      74             :         constexpr data_t blob_projected(data_t s, SelfType_t<data_t> a, SelfType_t<data_t> alpha,
      75             :                                         int m)
      76       94772 :         {
      77             :             // expect alpha > 0
      78       94772 :             using namespace elsa::math;
      79             : 
      80       94772 :             const data_t sda = s / a;
      81       94772 :             const data_t sdas = std::pow(sda, 2);
      82       94772 :             const data_t w = 1.0 - sdas;
      83             : 
      84       94772 :             if (w > 1.0e-10) {
      85       94460 :                 const auto arg = alpha * std::sqrt(w);
      86       94460 :                 if (m == 0) {
      87        2320 :                     return (2.0 * a / alpha) * std::sinh(arg) / math::bessi0(alpha);
      88       92140 :                 } else if (m == 1) {
      89        2320 :                     return (2.0 * a / alpha) * std::sqrt(w)
      90        2320 :                            * (std::cosh(arg) - std::sinh(arg) / arg) / math::bessi1(alpha);
      91             : 
      92       89820 :                 } else if (m == 2) {
      93       89820 :                     return (2.0 * a / alpha) * w
      94       89820 :                            * ((3.0 / (arg * arg) + 1.0) * std::sinh(arg) - (3.0 / arg) * cosh(arg))
      95       89820 :                            / bessi2(alpha);
      96       89820 :                 } else {
      97           0 :                     throw Error("m out of range in blob_projected()");
      98           0 :                 }
      99         312 :             }
     100         312 :             return 0.0;
     101         312 :         }
     102             : 
     103             :         template <typename data_t>
     104             :         constexpr data_t blob_projected(data_t s)
     105           0 :         {
     106           0 :             return blob_projected(s, 2.f, 10.83f, 2);
     107           0 :         }
     108             :     } // namespace blobs
     109             : 
     110             :     template <typename data_t>
     111             :     class Blob
     112             :     {
     113             :     public:
     114             :         constexpr Blob(data_t radius, SelfType_t<data_t> alpha, int order);
     115             : 
     116             :         constexpr data_t operator()(data_t s);
     117             : 
     118             :         constexpr data_t radius() const;
     119             : 
     120             :         constexpr data_t alpha() const;
     121             : 
     122             :         constexpr int order() const;
     123             : 
     124             :     private:
     125             :         data_t radius_;
     126             :         data_t alpha_;
     127             :         int order_;
     128             :     };
     129             : 
     130             :     template <typename data_t>
     131             :     class ProjectedBlob
     132             :     {
     133             :     public:
     134             :         constexpr ProjectedBlob(data_t radius, SelfType_t<data_t> alpha, int order)
     135             :             : radius_(radius), alpha_(alpha), order_(order)
     136         947 :         {
     137         947 :         }
     138             : 
     139             :         constexpr data_t operator()(data_t s)
     140       94772 :         {
     141       94772 :             return blobs::blob_projected(s, radius_, alpha_, order_);
     142       94772 :         }
     143             : 
     144       54062 :         constexpr data_t radius() const { return radius_; }
     145             : 
     146             :         constexpr data_t alpha() const { return alpha_; }
     147             : 
     148             :         constexpr int order() const { return order_; }
     149             : 
     150             :     private:
     151             :         data_t radius_;
     152             :         data_t alpha_;
     153             :         int order_;
     154             :     };
     155             : 
     156             : } // namespace elsa