Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include <optional>
       4             : 
       5             : #include "DataContainer.h"
       6             : #include "Solver.h"
       7             : #include "Scaling.h"
       8             : #include "MaybeUninitialized.hpp"
       9             : 
      10             : namespace elsa
      11             : {
      12             :     /// @brief Default constraint for Landweber Iterations, i.e. no constraint
      13             :     struct IdentityProjection {
      14             :         template <class data_t>
      15             :         void operator()([[maybe_unused]] DataContainer<data_t>& x) const
      16           4 :         {
      17           4 :         }
      18             :     };
      19             : 
      20             :     /// @brief Constraint functor, which constraints a `DataContainer` such that,
      21             :     /// it for all values \f$x_i\f$ it holds: \f$0 < x_i < 1\f$
      22             :     struct BoxConstraint {
      23             :         template <class data_t>
      24             :         void operator()(DataContainer<data_t>& x) const
      25             :         {
      26             :             x = clip(x, static_cast<data_t>(0), static_cast<data_t>(1));
      27             :         }
      28             :     };
      29             : 
      30             :     /// @brief Constraint functor, which constraints a `DataContainer` such that,
      31             :     /// it for all values \f$x_i\f$ it holds: \f$0 < x_i\f$
      32             :     struct NonNegativeConstraint {
      33             :         template <class data_t>
      34             :         void operator()(DataContainer<data_t>& x) const
      35             :         {
      36             :             // TODO: This could be optimized, but for now this works
      37             :             auto maxVal = x.maxElement();
      38             :             x = clip(x, static_cast<data_t>(0), maxVal);
      39             :         }
      40             :     };
      41             : 
      42             :     /**
      43             :      * @brief Base class for Landweber iteration. These are a group of algorithms solving ill-posed
      44             :      * linear inverse problem.
      45             :      *
      46             :      * The algorithms solve the equation of the form \f$A x = b\f$, by iteratively solving
      47             :      * \f$ \argmin_x \frac{1}{2} \| Ax - b \|_{2}^2 \f$. The general update rule for Landweber
      48             :      * iterations is:
      49             :      *
      50             :      * - \f$ x_{k+1} =  x_{k} + \lambda T A^T M (A(x_{k}) - b)\f$
      51             :      *
      52             :      * Here, \f$k\f$ is the number of the current iteration, \f$b\f$ the measured data (i.e.
      53             :      * sinogram), \f$A\f$ an operator (i.e. the Radon transform). \f$T\f$ and \f$M\f$ are method
      54             :      * specific, the exact choice of these determine the exact method. Among others:
      55             :      *
      56             :      * - Given \f$ T = M = I \f$ then the algorithm is referred to as the classical Landweber
      57             :      *   iteration
      58             :      * - Given \f$ T = I\f$ \f$ M = \frac{1}{m} \text{diag}(\frac{1}{ \| a_i \|^2_2 }) \f$ (where
      59             :      *   \f$ \| a_i \|^2_2 \f$ is the \f$\ell^2\f$-norm of the \f$i\f$-th row of the operator)
      60             :      *   referred to as Cimmino's method
      61             :      * - Given \f$ T = \text{diag}(\text{row sum})^{-1}\f$ \f$ M = \text{diag}(\text{col sum})^1\f$
      62             :      *   it is referred to as SART
      63             :      *
      64             :      * @author David Frank
      65             :      */
      66             :     template <typename data_t = real_t>
      67             :     class LandweberIteration : public Solver<data_t>
      68             :     {
      69             :     public:
      70             :         /// Scalar alias
      71             :         using Scalar = typename Solver<data_t>::Scalar;
      72             : 
      73             :         /**
      74             :          * @brief Constructor for Landweber type solver, accepting an operator, a measurement vector
      75             :          * and a step size#
      76             :          *
      77             :          * @param[in] A linear operator to solve the problem with
      78             :          * @param[in] b measurement vector of the problem
      79             :          * @param[in] stepSize the fixed step size to be used while solving
      80             :          */
      81             :         LandweberIteration(const LinearOperator<data_t>& A, const DataContainer<data_t>& b,
      82             :                            data_t stepSize);
      83             : 
      84             :         /**
      85             :          * @brief Constructor for Landweber type solver, accepting an operator and a measurement
      86             :          * vector
      87             :          *
      88             :          * @param[in] A linear operator to solve the problem with
      89             :          * @param[in] b measurement vector of the problem
      90             :          */
      91             :         LandweberIteration(const LinearOperator<data_t>& A, const DataContainer<data_t>& b);
      92             : 
      93             :         /// make copy constructor deletion explicit
      94             :         LandweberIteration(const LandweberIteration<data_t>&) = delete;
      95             : 
      96             :         /// default destructor
      97           8 :         ~LandweberIteration() override = default;
      98             : 
      99             :         void setProjection(const std::function<void(DataContainer<data_t>&)> projection);
     100             : 
     101             :         DataContainer<data_t> setup(std::optional<DataContainer<data_t>> x) override;
     102             : 
     103             :         DataContainer<data_t> step(DataContainer<data_t> x) override;
     104             : 
     105             :         std::string formatHeader() const override;
     106             : 
     107             :         std::string formatStep(const DataContainer<data_t>& x) const override;
     108             : 
     109             :     protected:
     110             :         /// Setup the \f$T * A^T * M\f& operator, implemented by the base classes to allow for
     111             :         /// different types of solvers
     112             :         virtual std::unique_ptr<LinearOperator<data_t>>
     113             :             setupOperators(const LinearOperator<data_t>& A) const = 0;
     114             : 
     115             :         /// the linear operator
     116             :         std::unique_ptr<LinearOperator<data_t>> A_;
     117             : 
     118             :         /// the measurement data
     119             :         DataContainer<data_t> b_;
     120             : 
     121             :         /// temporary for storage efficient computation of the residual
     122             :         DataContainer<data_t> residual_;
     123             : 
     124             :         std::function<void(DataContainer<data_t>&)> projection_ = IdentityProjection{};
     125             : 
     126             :         /// The composition of T, A^T and M to apply in each update
     127             :         std::unique_ptr<LinearOperator<data_t>> tam_{};
     128             : 
     129             :         /// the step size
     130             :         MaybeUninitialized<data_t> stepSize_;
     131             : 
     132             :         /// implement the polymorphic comparison operation
     133             :         bool isEqual(const Solver<data_t>& other) const override;
     134             :     };
     135             : } // namespace elsa