Line data    Source code

       1             : #include "AXDTOperator.h"
       2             : #include "SphericalCoefficientsDescriptor.h"
       3             : #include "WeightingFunction.h"
       4             : #include "BlockLinearOperator.h"
       5             : #include "DataContainer.h"
       6             : #include "Scaling.h"
       7             : #include "Timer.h"
       8             : #include "TypeCasts.hpp"
       9             : #include "elsaDefines.h"
      10             : 
      11             : namespace elsa
      12             : {
      13             : 
      14             :     template <typename data_t>
      15             :     AXDTOperator<data_t>::AXDTOperator(const SphericalCoefficientsDescriptor& domainDescriptor,
      16             :                                        const XGIDetectorDescriptor& rangeDescriptor,
      17             :                                        const LinearOperator<data_t>& projector,
      18             :                                        std::optional<DataContainer<data_t>> weights)
      19             :         : LinearOperator<data_t>(domainDescriptor, rangeDescriptor)
      20          12 :     {
      21             : 
      22          12 :         auto W = [&]() {
      23          12 :             if (weights)
      24           4 :                 return *weights;
      25           8 :             else
      26           8 :                 return axdt::exactWeightingFunction<data_t>(
      27           8 :                     rangeDescriptor, domainDescriptor.symmetry, domainDescriptor.degree);
      28          12 :         }();
      29             : 
      30             :         // Number of spherical harmonics
      31          12 :         assert(W.getNumberOfBlocks() == domainDescriptor.getNumberOfBlocks());
      32             : 
      33          12 :         OperatorList ops;
      34             : 
      35             :         // create composite operators of projector and scalings
      36         232 :         for (index_t i = 0; i < W.getNumberOfBlocks(); ++i) {
      37             :             // create a scaling operator
      38         220 :             Scaling<data_t> W_i{materialize(W.getBlock(i))};
      39             : 
      40         220 :             ops.emplace_back((W_i * projector).clone());
      41         220 :         }
      42             : 
      43          12 :         bl_op = std::make_unique<BlockLinearOperator<data_t>>(
      44          12 :             *this->_domainDescriptor, rangeDescriptor, ops,
      45          12 :             BlockLinearOperator<data_t>::BlockType::COL);
      46          12 :     }
      47             : 
      48             :     template <typename data_t>
      49             :     AXDTOperator<data_t>::AXDTOperator(const AXDTOperator& other)
      50             :         : LinearOperator<data_t>(*other._domainDescriptor, *other._rangeDescriptor)
      51           0 :     {
      52           0 :         bl_op = downcast<BlockLinearOperator<data_t>>(other.bl_op->clone());
      53           0 :     }
      54             : 
      55             :     template <typename data_t>
      56             :     AXDTOperator<data_t>* AXDTOperator<data_t>::cloneImpl() const
      57           0 :     {
      58           0 :         return new AXDTOperator<data_t>(*this);
      59           0 :     }
      60             : 
      61             :     template <typename data_t>
      62             :     bool AXDTOperator<data_t>::isEqual(const LinearOperator<data_t>& other) const
      63           0 :     {
      64           0 :         if (!LinearOperator<data_t>::isEqual(other))
      65           0 :             return false;
      66             : 
      67             :         // static_cast as type checked in base comparison
      68           0 :         const auto& otherOp = downcast<const AXDTOperator<data_t>>(other);
      69             : 
      70           0 :         return *bl_op == *(otherOp.bl_op);
      71           0 :     }
      72             : 
      73             :     template <typename data_t>
      74             :     void AXDTOperator<data_t>::applyImpl(const DataContainer<data_t>& x,
      75             :                                          DataContainer<data_t>& Ax) const
      76           8 :     {
      77           8 :         Timer timeguard("AXDTOperator", "apply");
      78             : 
      79           8 :         bl_op->apply(x, Ax);
      80           8 :     }
      81             : 
      82             :     template <typename data_t>
      83             :     void AXDTOperator<data_t>::applyAdjointImpl(const DataContainer<data_t>& y,
      84             :                                                 DataContainer<data_t>& Aty) const
      85           4 :     {
      86           4 :         Timer timeguard("AXDTOperator", "applyAdjoint");
      87             : 
      88           4 :         bl_op->applyAdjoint(y, Aty);
      89             :         //        Logger::get("AXDTOperator")->info("ApplyAdjoint result {}", Aty.sum());
      90           4 :     }
      91             : 
      92             :     template class AXDTOperator<float>;
      93             :     template class AXDTOperator<double>;
      94             : 
      95             : } // namespace elsa