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             :         constexpr inline index_t factorial(index_t n) noexcept
      11       43055 :         {
      12       43055 :             return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
      13       43055 :         }
      14             : 
      15             :         /// Compute binomial coefficient
      16             :         constexpr inline index_t binom(index_t n, index_t k) noexcept
      17       23057 :         {
      18       23057 :             return (k > n)
      19       23057 :                        ? 0
      20       23057 :                        : (k == 0 || k == n) ? 1
      21       23056 :                                             : (k == 1 || k == n - 1)
      22       10802 :                                                   ? n
      23       10802 :                                                   : (k + k < n) ? (binom(n - 1, k - 1) * n) / k
      24        1532 :                                                                 : (binom(n - 1, k) * n) / (n - k);
      25       23057 :         }
      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             :         constexpr data_t heaviside(data_t x1, data_t c)
      38       21704 :         {
      39       21704 :             if (x1 == 0) {
      40        1052 :                 return c;
      41       20652 :             } else if (x1 < 0) {
      42        9588 :                 return 0;
      43       11064 :             } else {
      44       11064 :                 return 1;
      45       11064 :             }
      46       21704 :         }
      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             :     data_t meyerFunction(data_t x)
      53      294762 :     {
      54      294762 :         if (x < 0) {
      55      156060 :             return 0;
      56      156060 :         } else if (0 <= x && x <= 1) {
      57      138702 :             return 35 * std::pow(x, 4) - 84 * std::pow(x, 5) + 70 * std::pow(x, 6)
      58      138702 :                    - 20 * std::pow(x, 7);
      59      138702 :         } else {
      60           0 :             return 1;
      61           0 :         }
      62      294762 :     }
      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             :         data_t b(data_t w)
      70      417384 :         {
      71      417384 :             if (1 <= std::abs(w) && std::abs(w) <= 2) {
      72       56484 :                 return std::sin(pi<data_t> / 2.0 * meyerFunction(std::abs(w) - 1));
      73      360900 :             } else if (2 < std::abs(w) && std::abs(w) <= 4) {
      74       29586 :                 return std::cos(pi<data_t> / 2.0 * meyerFunction(1.0 / 2 * std::abs(w) - 1));
      75      331314 :             } else {
      76      331314 :                 return 0;
      77      331314 :             }
      78      417384 :         }
      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             :         data_t phi(data_t w)
      84        6144 :         {
      85        6144 :             if (std::abs(w) <= 1.0 / 2) {
      86           6 :                 return 1;
      87        6138 :             } 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        6138 :             } else {
      90        6138 :                 return 0;
      91        6138 :             }
      92        6144 :         }
      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             :         data_t phiHat(data_t w, data_t h)
      98        6144 :         {
      99        6144 :             if (std::abs(h) <= std::abs(w)) {
     100        3258 :                 return phi(w);
     101        3258 :             } else {
     102        2886 :                 return phi(h);
     103        2886 :             }
     104        6144 :         }
     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             :         data_t psiHat1(data_t w)
     110      208692 :         {
     111      208692 :             return std::sqrt(std::pow(b(2 * w), 2) + std::pow(b(w), 2));
     112      208692 :         }
     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             :         data_t psiHat2(data_t w)
     118      208692 :         {
     119      208692 :             if (w <= 0) {
     120      106752 :                 return std::sqrt(meyerFunction(1 + w));
     121      106752 :             } else {
     122      101940 :                 return std::sqrt(meyerFunction(1 - w));
     123      101940 :             }
     124      208692 :         }
     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             :         data_t psiHat(data_t w, data_t h)
     130      208896 :         {
     131      208896 :             if (w == 0) {
     132         204 :                 return 0;
     133      208692 :             } else {
     134      208692 :                 return psiHat1(w) * psiHat2(h / w);
     135      208692 :             }
     136      208896 :         }
     137             :     } // namespace shearlet
     138             : } // namespace elsa