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