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