Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "LinearOperator.h"
       4             : 
       5             : #include <vector>
       6             : 
       7             : namespace elsa
       8             : {
       9             :     /**
      10             :      * @brief Operator to compute finite differences.
      11             :      *
      12             :      * @author Matthias Wieczorek - initial code
      13             :      * @author Maximilian Hornung - rewrite and performance optimization
      14             :      * @author Tobias Lasser - modernization
      15             :      *
      16             :      * @tparam data_t data type for the domain and range of the operator, defaulting to real_t
      17             :      *
      18             :      * This class represents a linear operator D that computes finite differences,
      19             :      * using the central, forward, or backward differences.
      20             :      */
      21             :     template <typename data_t = real_t>
      22             :     class FiniteDifferences : public LinearOperator<data_t>
      23             :     {
      24             :     public:
      25             :         /// supported types of finite differences
      26             :         enum class DiffType { FORWARD, BACKWARD, CENTRAL };
      27             : 
      28             :         /**
      29             :          * @brief Constructor for FiniteDifferences over all dimensions.
      30             :          *
      31             :          * @param[in] domainDescriptor Descriptor for domain
      32             :          * @param[in] type denoting the type of finite differences
      33             :          *
      34             :          * This implementation uses zero padding such that it's equal to the
      35             :          * following matrix formulation (in 1D)
      36             :          * - Dforward  = full(spdiags([-e e], 0:1, n, n));
      37             :          * - Dbackward = full(spdiags([-e e], -1:0, n, n));
      38             :          * - Dcentral  = spdiags(0.5*[-e e], [-1,1], n, n);
      39             :          *
      40             :          * Note: the descriptor for the range is automatically generated from the domain.
      41             :          */
      42             :         explicit FiniteDifferences(const DataDescriptor& domainDescriptor,
      43             :                                    DiffType type = DiffType::FORWARD);
      44             : 
      45             :         /**
      46             :          * @brief Constructor for FiniteDifferences over selected dimensions.
      47             :          *
      48             :          * @param[in] domainDescriptor Descriptor for domain
      49             :          * @param[in] activeDims Boolean vector defining which dimensions are active or not
      50             :          * @param[in] type denoting the type of finite differences
      51             :          *
      52             :          * This implementation uses zero padding such that it's equal to the
      53             :          * following matrix formulation (in 1D)
      54             :          * - Dforward  = full(spdiags([-e e], 0:1, n, n));
      55             :          * - Dbackward = full(spdiags([-e e], -1:0, n, n));
      56             :          * - Dcentral  = spdiags(0.5*[-e e], [-1,1], n, n);
      57             :          *
      58             :          * Note: the descriptor for the range is automatically generated from the domain.
      59             :          */
      60             :         FiniteDifferences(const DataDescriptor& domainDescriptor, const BooleanVector_t& activeDims,
      61             :                           DiffType type = DiffType::FORWARD);
      62             : 
      63             :         /// default destructor
      64           0 :         ~FiniteDifferences() override = default;
      65             : 
      66             :     protected:
      67             :         /// default copy constructor, hidden from non-derived classes to prevent potential slicing
      68             :         FiniteDifferences(const FiniteDifferences<data_t>&) = default;
      69             : 
      70             :         /// apply the finite differences operator
      71             :         void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
      72             : 
      73             :         /// apply the adjoint of the finite differences operator
      74             :         void applyAdjointImpl(const DataContainer<data_t>& y,
      75             :                               DataContainer<data_t>& Aty) const override;
      76             : 
      77             :         /// implement the polymorphic clone operation
      78             :         FiniteDifferences<data_t>* cloneImpl() const override;
      79             : 
      80             :         /// implement the polymorphic comparison operation
      81             :         bool isEqual(const LinearOperator<data_t>& other) const override;
      82             : 
      83             :     private:
      84             :         /// type of the finite differences to be computed
      85             :         DiffType _type;
      86             : 
      87             :         /// boolean vector for active dimensions when computing finite differences
      88             :         BooleanVector_t _activeDims;
      89             : 
      90             :         /// precompute some helper variables to optimize speed
      91             :         void precomputeHelpers();
      92             : 
      93             :         IndexVector_t _coordDiff;  /// precomputed helper for coordinate diffs
      94             :         IndexVector_t _coordDelta; /// precomputed helper for coordinate deltas
      95             :         IndexVector_t _dimCounter; /// precomputed helper for active dim counter
      96             : 
      97             :         /// the actual finite differences computations (with mode as template parameter for
      98             :         /// performance)
      99             :         template <typename FDtype>
     100             :         void applyHelper(const DataContainer<data_t>& x, DataContainer<data_t>& Ax,
     101             :                          FDtype type) const;
     102             : 
     103             :         /// the actual finite differences computations (with mode as template parameter for
     104             :         /// performance)
     105             :         template <typename FDtype>
     106             :         void applyAdjointHelper(const DataContainer<data_t>& y, DataContainer<data_t>& Aty,
     107             :                                 FDtype type) const;
     108             :     };
     109             : } // namespace elsa