Line data    Source code

       1             : #include "SiddonsMethod.h"
       2             : #include "Timer.h"
       3             : #include "TraverseAABB.h"
       4             : #include "TypeCasts.hpp"
       5             : 
       6             : #include <stdexcept>
       7             : #include <type_traits>
       8             : 
       9             : namespace elsa
      10             : {
      11             :     template <typename data_t>
      12           0 :     SiddonsMethod<data_t>::SiddonsMethod(const VolumeDescriptor& domainDescriptor,
      13             :                                          const DetectorDescriptor& rangeDescriptor)
      14           0 :         : base_type(domainDescriptor, rangeDescriptor)
      15             :     {
      16           0 :         auto dim = domainDescriptor.getNumberOfDimensions();
      17           0 :         if (dim != rangeDescriptor.getNumberOfDimensions()) {
      18           0 :             throw LogicError("SiddonsMethod: domain and range dimension need to match");
      19             :         }
      20             : 
      21           0 :         if (dim != 2 && dim != 3) {
      22           0 :             throw LogicError("SiddonsMethod: only supporting 2d/3d operations");
      23             :         }
      24             : 
      25           0 :         if (rangeDescriptor.getNumberOfGeometryPoses() == 0) {
      26           0 :             throw LogicError("SiddonsMethod: geometry list was empty");
      27             :         }
      28           0 :     }
      29             : 
      30             :     template <typename data_t>
      31           0 :     SiddonsMethod<data_t>* SiddonsMethod<data_t>::_cloneImpl() const
      32             :     {
      33           0 :         return new self_type(downcast<VolumeDescriptor>(*this->_domainDescriptor),
      34           0 :                              downcast<DetectorDescriptor>(*this->_rangeDescriptor));
      35             :     }
      36             : 
      37             :     template <typename data_t>
      38           0 :     bool SiddonsMethod<data_t>::_isEqual(const LinearOperator<data_t>& other) const
      39             :     {
      40           0 :         if (!LinearOperator<data_t>::isEqual(other))
      41           0 :             return false;
      42             : 
      43           0 :         auto otherSM = downcast_safe<SiddonsMethod>(&other);
      44           0 :         return static_cast<bool>(otherSM);
      45             :     }
      46             : 
      47             :     template <typename data_t>
      48           0 :     data_t SiddonsMethod<data_t>::traverseRayForward(BoundingBox aabb, const RealRay_t& ray,
      49             :                                                      const DataContainer<data_t>& x) const
      50             :     {
      51           0 :         const auto& domain = x.getDataDescriptor();
      52             : 
      53             :         // --> setup traversal algorithm
      54           0 :         TraverseAABB traverse(aabb, ray);
      55             : 
      56           0 :         data_t accumulator = data_t(0);
      57             : 
      58             :         // --> initial index to access the data vector
      59           0 :         auto dataIndexForCurrentVoxel = domain.getIndexFromCoordinate(traverse.getCurrentVoxel());
      60             : 
      61           0 :         while (traverse.isInBoundingBox()) {
      62             : 
      63           0 :             auto weight = traverse.updateTraverseAndGetDistance();
      64             :             // --> update result depending on the operation performed
      65           0 :             accumulator += x[dataIndexForCurrentVoxel] * weight;
      66             : 
      67           0 :             dataIndexForCurrentVoxel = domain.getIndexFromCoordinate(traverse.getCurrentVoxel());
      68             :         }
      69             : 
      70           0 :         return accumulator;
      71           0 :     }
      72             : 
      73             :     template <typename data_t>
      74           0 :     void SiddonsMethod<data_t>::traverseRayBackward(BoundingBox aabb, const RealRay_t& ray,
      75             :                                                     const value_type& detectorValue,
      76             :                                                     DataContainer<data_t>& Aty) const
      77             :     {
      78           0 :         const auto& domain = Aty.getDataDescriptor();
      79             : 
      80             :         // --> setup traversal algorithm
      81           0 :         TraverseAABB traverse(aabb, ray);
      82             : 
      83             :         // --> initial index to access the data vector
      84           0 :         auto dataIndexForCurrentVoxel = domain.getIndexFromCoordinate(traverse.getCurrentVoxel());
      85             : 
      86           0 :         while (traverse.isInBoundingBox()) {
      87             : 
      88           0 :             auto weight = traverse.updateTraverseAndGetDistance();
      89             : 
      90           0 : #pragma omp atomic
      91           0 :             Aty[dataIndexForCurrentVoxel] += detectorValue * weight;
      92             : 
      93           0 :             dataIndexForCurrentVoxel = domain.getIndexFromCoordinate(traverse.getCurrentVoxel());
      94             :         }
      95           0 :     }
      96             : 
      97             :     // ------------------------------------------
      98             :     // explicit template instantiation
      99             :     template class SiddonsMethod<float>;
     100             :     template class SiddonsMethod<double>;
     101             : } // namespace elsa