Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "AXDTOperator.h"
       4             : #include "Functional.h"
       5             : #include "RandomBlocksDescriptor.h"
       6             : 
       7             : namespace elsa
       8             : {
       9             :     /**
      10             :      * @brief Implement the loss functions presented by Schilling et al (see ref below),
      11             :      * which implement the log-likelihood functions for AXDT reconstruction, which
      12             :      * follow the Rician distribution.
      13             :      *
      14             :      * We assume the measurements follow the following formula:
      15             :      * \[
      16             :      * s_{n,j} = a_j + b_j * \cos(\frac{2 \pi n}{N} - \phi_j),
      17             :      * \]
      18             :      * This function assumes that \f$ A_j \sim \mathcal{N}(a_j, \frac{a_j}{N})\$, and \f$
      19             :      * \frac{1}{2} B_j \sim \mathcal{R} \left( \frac{b_j}{2}, \sqrt{\frac{a_j}{2N}} \right) \f$.
      20             :      * Here, capital letters represent the random variables. This formula is an
      21             :      * accurate assumption for the measurement of AXDT.
      22             :      *
      23             :      * Schilling derived the proper log likelihood functionals for the above assumptions.
      24             :      * They also present two functionals, one where the Rician is approximated
      25             :      * by a Gaussian, which holds under certain assumptions. For the exact formulation,
      26             :      * derivative and Hessian, please check the references below.
      27             :      *
      28             :      * References:
      29             :      * - Statistical Models for Anisotropic X-Ray Dark-field Tomography, Schilling et al (2017)
      30             :      *
      31             :      */
      32             :     template <typename data_t = real_t>
      33             :     class RicianLoss : public Functional<data_t>
      34             :     {
      35             :     public:
      36             :         /**
      37             :          * @brief Construct a Joint Rician Loss
      38             :          * the domain consists of the cartesian product of the attenuation data and the spherical
      39             :          * harmonics coefficients
      40             :          *
      41             :          * @param[in] ffa the flat-field measured value of RV a
      42             :          * @param[in] ffb the flat-field measured value of RV b
      43             :          * @param[in] a the measured value of RV a (with the measured sample in-place)
      44             :          * @param[in] b the measured value of RV b (with the measured sample in-place)
      45             :          * @param[in] absorp_op operator modeling the projection of the regular absorption signal
      46             :          * @param[in] axdt_op operator modeling the projection of the dark-field signal
      47             :          * @param[in] N total grating-stepping steps
      48             :          * @param[in] approximate if true, use a Gaussian distribution to approximate the Racian
      49             :          */
      50             :         RicianLoss(const DataContainer<data_t>& ffa, const DataContainer<data_t>& ffb,
      51             :                    const DataContainer<data_t>& a, const DataContainer<data_t>& b,
      52             :                    const LinearOperator<data_t>& absorp_op, const LinearOperator<data_t>& axdt_op,
      53             :                    index_t N, bool approximate);
      54             : 
      55             :         /// make copy constructor deletion explicit
      56             :         RicianLoss(const RicianLoss<data_t>&) = delete;
      57             : 
      58             :         /// default destructor
      59          24 :         ~RicianLoss() override = default;
      60             : 
      61             :         bool isDifferentiable() const override;
      62             : 
      63             :     protected:
      64             :         /// the evaluation of this functional
      65             :         data_t evaluateImpl(const DataContainer<data_t>& Rx) const override;
      66             : 
      67             :         /// the computation of the gradient
      68             :         void getGradientImpl(const DataContainer<data_t>& Rx,
      69             :                              DataContainer<data_t>& out) const override;
      70             : 
      71             :         /// the computation of the Hessian
      72             :         LinearOperator<data_t> getHessianImpl(const DataContainer<data_t>& Rx) const override;
      73             : 
      74             :         /// implement the polymorphic clone operation
      75             :         RicianLoss<data_t>* cloneImpl() const override;
      76             : 
      77             :         /// implement the polymorphic comparison operation
      78             :         bool isEqual(const Functional<data_t>& other) const override;
      79             : 
      80             :     private:
      81             :         /// the flat-field measured value of RV a
      82             :         DataContainer<data_t> ffa_;
      83             : 
      84             :         /// the flat-field measured value of RV b
      85             :         DataContainer<data_t> ffb_;
      86             : 
      87             :         /// the measured value of RV a
      88             :         DataContainer<data_t> a_tilde_;
      89             : 
      90             :         /// the measured value of RV b
      91             :         DataContainer<data_t> b_tilde_;
      92             : 
      93             :         /// operator modeling the projection of the regular absorption signal
      94             :         std::unique_ptr<LinearOperator<data_t>> absorp_op_;
      95             : 
      96             :         /// operator modeling the projection of the dark-field signal
      97             :         std::unique_ptr<LinearOperator<data_t>> axdt_op_;
      98             : 
      99             :         /// store the frequently used value alpha (= ffb / ffa)
     100             :         DataContainer<data_t> alpha_;
     101             : 
     102             :         /// cache the logarithm to save unnecessary computations
     103             :         DataContainer<data_t> log_ffa_;
     104             : 
     105             :         /// total grating-stepping steps
     106             :         data_t N_;
     107             : 
     108             :         /// bool to indicate whether to use a Gaussian distribution to approximate the Racian
     109             :         bool approximate_;
     110             :     };
     111             : 
     112             : } // namespace elsa