Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "DataContainer.h"
       4             : #include "DataDescriptor.h"
       5             : #include "LinearOperator.h"
       6             : #include "PartitionDescriptor.h"
       7             : #include "Scaling.h"
       8             : #include "TypeTraits.hpp"
       9             : #include "VolumeDescriptor.h"
      10             : #include "elsaDefines.h"
      11             : #include <cstdlib>
      12             : #include <functional>
      13             : 
      14             : namespace elsa
      15             : {
      16             : 
      17             :     template <typename data_t = float>
      18             :     using Filter = Scaling<add_complex_t<data_t>>;
      19             : 
      20             :     /**
      21             :      * @brief Generates filter form transfer function in k-Space
      22             :      *
      23             :      * @tparam data_t underlying floating point type
      24             :      * @param descriptor Data descriptor describing filter dimensions
      25             :      * @param transferFunction Transfer function lambda receiving a wavenumber vector
      26             :      * @return Filter<data_t> Scaling with DC bin at position 0
      27             :      */
      28             :     template <typename data_t = float>
      29             :     std::unique_ptr<Filter<data_t>>
      30             :         makeFilter(const DataDescriptor& descriptor,
      31             :                    std::function<add_complex_t<data_t>(IndexVector_t)> transferFunction)
      32           5 :     {
      33             : 
      34           5 :         auto dim = descriptor.getNumberOfDimensions();
      35           5 :         auto coefficients = descriptor.getNumberOfCoefficientsPerDimension();
      36           5 :         coefficients[dim - 1] = 1;
      37             : 
      38           5 :         auto sliceDesc = VolumeDescriptor{coefficients};
      39           5 :         auto midPoint = coefficients / 2;
      40             : 
      41           5 :         DataContainer<add_complex_t<data_t>> filter{sliceDesc};
      42             : 
      43           5 :         if (!(dim == 2 || dim == 3)) {
      44           0 :             throw InvalidArgumentError("makeFilter:: Should only be used for 2D or 3D sinograms");
      45           0 :         }
      46             : 
      47           5 : #pragma omp parallel for
      48           5 :         for (index_t i = 0; i < sliceDesc.getNumberOfCoefficients(); ++i) {
      49           0 :             IndexVector_t nu = sliceDesc.getCoordinateFromIndex(i) - midPoint;
      50           0 :             auto H = transferFunction(nu);
      51           0 :             filter[i] = H;
      52           0 :         }
      53             : 
      54           5 :         return std::make_unique<Filter<data_t>>(ifftShift(filter));
      55           5 :     }
      56             : 
      57             :     template <typename data_t = float>
      58             :     data_t kMax(const DataDescriptor& descriptor)
      59           5 :     {
      60           5 :         return 2 * pi_t
      61           5 :                * (descriptor.getNumberOfCoefficientsPerDimension().head(
      62           5 :                       descriptor.getNumberOfDimensions() - 1)
      63           5 :                   / 2)
      64           5 :                      .norm();
      65           5 :     }
      66             : 
      67             :     // TODO: Non-mononotonic filter function start to increase again for 2d filters (k_max  *=
      68             :     // sqrt(2)?)
      69             :     template <typename data_t = float>
      70             :     std::unique_ptr<Filter<data_t>> makeRamLakFilter(const DataDescriptor& descriptor)
      71           5 :     {
      72           5 :         auto n = descriptor.getNumberOfCoefficientsPerDimension()[0] / 2;
      73           5 :         auto kmax = kMax<data_t>(descriptor);
      74          81 :         return makeFilter<data_t>(descriptor, [&](IndexVector_t nu) {
      75          81 :             auto deltak = 2 * pi_t / n; // TODO: DC bin might have different factor in 3d
      76          81 :             auto k = 2 * pi_t * nu.norm();
      77          81 :             return (nu.isZero() ? 0.25 * deltak : k) / kmax;
      78          81 :         });
      79           5 :     }
      80             : 
      81             :     template <typename data_t = float>
      82             :     std::unique_ptr<Filter<data_t>> makeSheppLoganFilter(const DataDescriptor& descriptor)
      83             :     {
      84             :         auto kmax = kMax<data_t>(descriptor);
      85             :         return makeFilter<data_t>(descriptor, [&](IndexVector_t nu) {
      86             :             auto k = 2 * pi_t * nu.norm();
      87             :             return std::sin(pi_t / 2 * k / kmax) * 2 / pi_t;
      88             :         });
      89             :     }
      90             : 
      91             :     template <typename data_t = float>
      92             :     std::unique_ptr<Filter<data_t>> makeCosineFilter(const DataDescriptor& descriptor)
      93             :     {
      94             :         auto kmax = kMax<data_t>(descriptor);
      95             :         return makeFilter<data_t>(descriptor, [&](IndexVector_t nu) {
      96             :             auto k = 2 * pi_t * nu.norm();
      97             :             return std::sin(pi_t * k / kmax) / pi_t;
      98             :         });
      99             :     }
     100             : 
     101             :     template <typename data_t = float>
     102             :     std::unique_ptr<Filter<data_t>> makeHannFilter(const DataDescriptor& descriptor)
     103             :     {
     104             :         auto kmax = kMax<data_t>(descriptor);
     105             :         return makeFilter<data_t>(descriptor, [&](IndexVector_t nu) {
     106             :             auto k = 2 * pi_t * nu.norm();
     107             :             return k / kmax * std::cos(pi_t / 2 * k / kmax) * std::cos(pi_t / 2 * k / kmax);
     108             :         });
     109             :     }
     110             : 
     111             : } // namespace elsa