Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "Functional.h"
       4             : #include "RandomBlocksDescriptor.h"
       5             : 
       6             : namespace elsa
       7             : {
       8             :     /**
       9             :      * @brief Implement the loss functions presented by Schilling et al (see ref below),
      10             :      * which implement the log-likelihood functions for AXDT reconstruction, which
      11             :      * follow the Rician distribution.
      12             :      *
      13             :      * We assume the measurements follow the following formula:
      14             :      * \[
      15             :      * s_{n,j} = a_j + b_j * \cos(\frac{2 \pi n}{N} - \phi_j),
      16             :      * \]
      17             :      * This function assumes that \f$ A_j \sim \mathcal{N}(a_j, \frac{a_j}{N})\$, and \f$
      18             :      * \frac{1}{2} B_j \sim \mathcal{R} \left( \frac{b_j}{2}, \sqrt{\frac{a_j}{2N}} \right) \f$.
      19             :      * Here, capital letters represent the random variables. This formula is an
      20             :      * accurate assumption for the measurement of AXDT.
      21             :      *
      22             :      * Schilling derived the proper log likelihood functionals for the above assumptions.
      23             :      * They also present two functionals, one where the Rician is approximated
      24             :      * by a Gaussian, which holds under certain assumptions. For the exact formulation,
      25             :      * derivative and Hessian, please check the references below.
      26             :      *
      27             :      * References:
      28             :      * - Statistical Models for Anisotropic X-Ray Dark-field Tomography, Schilling et al (2017)
      29             :      *
      30             :      */
      31             :     template <typename data_t = real_t>
      32             :     class RicianLoss : public Functional<data_t>
      33             :     {
      34             :     public:
      35             :         /**
      36             :          * @brief full constructor for the AXDTStatRecon functional,
      37             :          * mapping domain vector to a scalar (without a residual)
      38             :          *
      39             :          * @param[in] ffa the flat-field measured value of RV a
      40             :          * @param[in] ffb the flat-field measured value of RV b
      41             :          * @param[in] a the measured value of RV a (with the measured sample in-place)
      42             :          * @param[in] b the measured value of RV b (with the measured sample in-place)
      43             :          * @param[in] absorp_op operator modeling the projection of the regular absorption signal
      44             :          * @param[in] axdt_op operator modeling the projection of the dark-field signal
      45             :          * @param[in] N total grating-stepping steps
      46             :          * @param[in] approximate if true, use a Gaussian distribution to approximate the Racian
      47             :          */
      48             :         RicianLoss(const DataContainer<data_t>& ffa, const DataContainer<data_t>& ffb,
      49             :                    const DataContainer<data_t>& a, const DataContainer<data_t>& b,
      50             :                    const LinearOperator<data_t>& absorp_op, const LinearOperator<data_t>& axdt_op,
      51             :                    index_t N, bool approximate);
      52             : 
      53             :         /// make copy constructor deletion explicit
      54             :         RicianLoss(const RicianLoss<data_t>&) = delete;
      55             : 
      56             :         /// default destructor
      57          24 :         ~RicianLoss() override = default;
      58             : 
      59             :         bool isDifferentiable() const override;
      60             : 
      61             :     protected:
      62             :         /// the evaluation of this functional
      63             :         data_t evaluateImpl(const DataContainer<data_t>& Rx) const override;
      64             : 
      65             :         /// the computation of the gradient
      66             :         void getGradientImpl(const DataContainer<data_t>& Rx,
      67             :                              DataContainer<data_t>& out) const override;
      68             : 
      69             :         /// the computation of the Hessian
      70             :         LinearOperator<data_t> getHessianImpl(const DataContainer<data_t>& Rx) const override;
      71             : 
      72             :         /// implement the polymorphic clone operation
      73             :         RicianLoss<data_t>* cloneImpl() const override;
      74             : 
      75             :         /// implement the polymorphic comparison operation
      76             :         bool isEqual(const Functional<data_t>& other) const override;
      77             : 
      78             :         /// generate a minimal volume descriptor as a placeholder when absorption data missing
      79             :         static std::unique_ptr<DataDescriptor> generate_placeholder_descriptor();
      80             : 
      81             :         /// construct a RandomBlockDescriptor to wrap up the absorption projection data
      82             :         /// and the axdt projection data as a single input for this functional
      83             :         static RandomBlocksDescriptor generate_descriptors(const DataDescriptor& desc1,
      84             :                                                            const DataDescriptor& desc2);
      85             : 
      86             :     private:
      87             :         /// the flat-field measured value of RV a
      88             :         DataContainer<data_t> ffa_;
      89             : 
      90             :         /// the flat-field measured value of RV b
      91             :         DataContainer<data_t> ffb_;
      92             : 
      93             :         /// the measured value of RV a
      94             :         DataContainer<data_t> a_tilde_;
      95             : 
      96             :         /// the measured value of RV b
      97             :         DataContainer<data_t> b_tilde_;
      98             : 
      99             :         /// operator modeling the projection of the regular absorption signal
     100             :         std::unique_ptr<LinearOperator<data_t>> absorp_op_;
     101             : 
     102             :         /// operator modeling the projection of the dark-field signal
     103             :         std::unique_ptr<LinearOperator<data_t>> axdt_op_;
     104             : 
     105             :         /// store the frequently used value alpha (= ffb / ffa)
     106             :         DataContainer<data_t> alpha_;
     107             : 
     108             :         /// store the frequently used value d_tilde (= b_tilde / a_tilde / alpha)
     109             :         DataContainer<data_t> d_tilde_;
     110             : 
     111             :         /// total grating-stepping steps
     112             :         data_t N_;
     113             : 
     114             :         /// bool to indicate whether to use a Gaussian distribution to approximate the Racian
     115             :         bool approximate_;
     116             :     };
     117             : 
     118             :     namespace axdt
     119             :     {
     120             :         /// element-wise logarithm of modified Bessel function of the first
     121             :         /// kind of order zero. i.e. log(B_i(0, x))
     122             :         template <typename data_t>
     123             :         DataContainer<data_t> log_bessel_0(const DataContainer<data_t>& x)
     124           4 :         {
     125           4 :             return bessel_log_0(x);
     126           4 :         }
     127             : 
     128             :         /// element-wise quotient between modified Bessel function of the first
     129             :         /// kind of order one and that of order zero. i.e. B_i(1, x) / B_i(0, x)
     130             :         template <typename data_t>
     131             :         DataContainer<data_t> quot_bessel_1_0(const DataContainer<data_t>& x)
     132           8 :         {
     133           8 :             return bessel_1_0(x);
     134           8 :         }
     135             : 
     136             :     } // namespace axdt
     137             : } // namespace elsa