Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "elsaDefines.h"
       4             : 
       5             : namespace elsa
       6             : {
       7             :     namespace math
       8             :     {
       9             :         /// Compute factorial \f$n!\f$ recursively
      10           0 :         constexpr inline index_t factorial(index_t n) noexcept
      11             :         {
      12           0 :             return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
      13             :         }
      14             : 
      15             :         /// Compute binomial coefficient
      16           0 :         constexpr inline index_t binom(index_t n, index_t k) noexcept
      17             :         {
      18             :             return (k > n)
      19           0 :                        ? 0
      20           0 :                        : (k == 0 || k == n) ? 1
      21           0 :                                             : (k == 1 || k == n - 1)
      22           0 :                                                   ? n
      23           0 :                                                   : (k + k < n) ? (binom(n - 1, k - 1) * n) / k
      24           0 :                                                                 : (binom(n - 1, k) * n) / (n - k);
      25             :         }
      26             : 
      27             :         /// Compute Heaviside-function
      28             :         /// \f[
      29             :         /// x \mapsto
      30             :         /// \begin{cases}
      31             :         /// 0: & x < 0 \\
      32             :         /// c: & x = 0 \\
      33             :         /// 1: & x > 0
      34             :         /// \end{cases}
      35             :         /// \f]
      36             :         template <typename data_t>
      37           0 :         constexpr data_t heaviside(data_t x1, data_t c)
      38             :         {
      39           0 :             if (x1 == 0) {
      40           0 :                 return c;
      41           0 :             } else if (x1 < 0) {
      42           0 :                 return 0;
      43             :             } else {
      44           0 :                 return 1;
      45             :             }
      46             :         }
      47             :     } // namespace math
      48             : 
      49             :     /// proposed in Y. Meyer, Oscillating Patterns in Image Processing and Nonlinear Evolution
      50             :     /// Equations. AMS, 2001
      51             :     template <typename data_t>
      52           0 :     data_t meyerFunction(data_t x)
      53             :     {
      54           0 :         if (x < 0) {
      55           0 :             return 0;
      56           0 :         } else if (0 <= x && x <= 1) {
      57           0 :             return 35 * std::pow(x, 4) - 84 * std::pow(x, 5) + 70 * std::pow(x, 6)
      58           0 :                    - 20 * std::pow(x, 7);
      59             :         } else {
      60           0 :             return 1;
      61             :         }
      62             :     }
      63             : 
      64             :     namespace shearlet
      65             :     {
      66             :         /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
      67             :         /// tutorial, 2014
      68             :         template <typename data_t>
      69           0 :         data_t b(data_t w)
      70             :         {
      71           0 :             if (1 <= std::abs(w) && std::abs(w) <= 2) {
      72           0 :                 return std::sin(pi<data_t> / 2.0 * meyerFunction(std::abs(w) - 1));
      73           0 :             } else if (2 < std::abs(w) && std::abs(w) <= 4) {
      74           0 :                 return std::cos(pi<data_t> / 2.0 * meyerFunction(1.0 / 2 * std::abs(w) - 1));
      75             :             } else {
      76           0 :                 return 0;
      77             :             }
      78             :         }
      79             : 
      80             :         /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
      81             :         /// tutorial, 2014
      82             :         template <typename data_t>
      83           0 :         data_t phi(data_t w)
      84             :         {
      85           0 :             if (std::abs(w) <= 1.0 / 2) {
      86           0 :                 return 1;
      87           0 :             } else if (1.0 / 2 < std::abs(w) && std::abs(w) < 1) {
      88           0 :                 return std::cos(pi<data_t> / 2.0 * meyerFunction(2 * std::abs(w) - 1));
      89             :             } else {
      90           0 :                 return 0;
      91             :             }
      92             :         }
      93             : 
      94             :         /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
      95             :         /// tutorial, 2014
      96             :         template <typename data_t>
      97           0 :         data_t phiHat(data_t w, data_t h)
      98             :         {
      99           0 :             if (std::abs(h) <= std::abs(w)) {
     100           0 :                 return phi(w);
     101             :             } else {
     102           0 :                 return phi(h);
     103             :             }
     104             :         }
     105             : 
     106             :         /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
     107             :         /// tutorial, 2014
     108             :         template <typename data_t>
     109           0 :         data_t psiHat1(data_t w)
     110             :         {
     111           0 :             return std::sqrt(std::pow(b(2 * w), 2) + std::pow(b(w), 2));
     112             :         }
     113             : 
     114             :         /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
     115             :         /// tutorial, 2014
     116             :         template <typename data_t>
     117           0 :         data_t psiHat2(data_t w)
     118             :         {
     119           0 :             if (w <= 0) {
     120           0 :                 return std::sqrt(meyerFunction(1 + w));
     121             :             } else {
     122           0 :                 return std::sqrt(meyerFunction(1 - w));
     123             :             }
     124             :         }
     125             : 
     126             :         /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
     127             :         /// tutorial, 2014
     128             :         template <typename data_t>
     129           0 :         data_t psiHat(data_t w, data_t h)
     130             :         {
     131           0 :             if (w == 0) {
     132           0 :                 return 0;
     133             :             } else {
     134           0 :                 return psiHat1(w) * psiHat2(h / w);
     135             :             }
     136             :         }
     137             :     } // namespace shearlet
     138             : } // namespace elsa