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          69 :     {
      16          69 :         auto dim = domainDescriptor.getNumberOfDimensions();
      17          69 :         if (dim != rangeDescriptor.getNumberOfDimensions()) {
      18           0 :             throw LogicError("SiddonsMethod: domain and range dimension need to match");
      19           0 :         }
      20             : 
      21          69 :         if (dim != 2 && dim != 3) {
      22           0 :             throw LogicError("SiddonsMethod: only supporting 2d/3d operations");
      23           0 :         }
      24             : 
      25          69 :         if (rangeDescriptor.getNumberOfGeometryPoses() == 0) {
      26           0 :             throw LogicError("SiddonsMethod: geometry list was empty");
      27           0 :         }
      28          69 :     }
      29             : 
      30             :     template <typename data_t>
      31             :     SiddonsMethod<data_t>* SiddonsMethod<data_t>::_cloneImpl() const
      32           2 :     {
      33           2 :         return new self_type(downcast<VolumeDescriptor>(*this->_domainDescriptor),
      34           2 :                              downcast<DetectorDescriptor>(*this->_rangeDescriptor));
      35           2 :     }
      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       33197 :     {
      51       33197 :         const auto& domain = x.getDataDescriptor();
      52             : 
      53       33197 :         if (domain.getNumberOfDimensions() == 2) {
      54       21639 :             return doTraverseRayForward<2>(aabb, ray, x, domain);
      55       21639 :         } else if (domain.getNumberOfDimensions() == 3) {
      56       11568 :             return doTraverseRayForward<3>(aabb, ray, x, domain);
      57       11568 :         }
      58             : 
      59 >1844*10^16 :         return data_t(0);
      60 >1844*10^16 :     }
      61             : 
      62             :     template <typename data_t>
      63             :     void SiddonsMethod<data_t>::traverseRayBackward(BoundingBox aabb, const RealRay_t& ray,
      64             :                                                     const value_type& detectorValue,
      65             :                                                     DataContainer<data_t>& Aty) const
      66       18988 :     {
      67       18988 :         const auto& domain = Aty.getDataDescriptor();
      68             : 
      69       18988 :         if (domain.getNumberOfDimensions() == 2) {
      70       18969 :             doTraverseRayBackward<2>(aabb, ray, detectorValue, Aty);
      71       18969 :         } else if (domain.getNumberOfDimensions() == 3) {
      72          14 :             doTraverseRayBackward<3>(aabb, ray, detectorValue, Aty);
      73          14 :         }
      74       18988 :     }
      75             : 
      76             :     template <typename data_t>
      77             :     template <int dim>
      78             :     data_t SiddonsMethod<data_t>::doTraverseRayForward(BoundingBox aabb, const RealRay_t& ray,
      79             :                                                        const DataContainer<data_t>& x,
      80             :                                                        const DataDescriptor& domain) const
      81       33191 :     {
      82             :         // --> setup traversal algorithm
      83       33191 :         TraverseAABB<dim> traverse(aabb, ray, domain.getProductOfCoefficientsPerDimension());
      84             : 
      85       33191 :         data_t accumulator = data_t(0);
      86             : 
      87             :         // --> initial index to access the data vector
      88       33191 :         auto dataIndexForCurrentVoxel = traverse.getCurrentIndex();
      89             : 
      90      970441 :         while (traverse.isInBoundingBox()) {
      91             : 
      92      937250 :             auto weight = traverse.updateTraverseAndGetDistance();
      93             :             // --> update result depending on the operation performed
      94      937250 :             accumulator += x[dataIndexForCurrentVoxel] * weight;
      95             : 
      96      937250 :             dataIndexForCurrentVoxel = traverse.getCurrentIndex();
      97      937250 :         }
      98             : 
      99       33191 :         return accumulator;
     100       33191 :     }
     101             : 
     102             :     template <typename data_t>
     103             :     template <int dim>
     104             :     void SiddonsMethod<data_t>::doTraverseRayBackward(BoundingBox aabb, const RealRay_t& ray,
     105             :                                                       const value_type& detectorValue,
     106             :                                                       DataContainer<data_t>& Aty) const
     107       19001 :     {
     108       19001 :         const auto& domain = Aty.getDataDescriptor();
     109             : 
     110             :         // --> setup traversal algorithm
     111       19001 :         TraverseAABB<dim> traverse(aabb, ray, domain.getProductOfCoefficientsPerDimension());
     112             :         // --> initial index to access the data vector
     113       19001 :         auto dataIndexForCurrentVoxel = traverse.getCurrentIndex();
     114             : 
     115      855789 :         while (traverse.isInBoundingBox()) {
     116             : 
     117      836788 :             auto weight = traverse.updateTraverseAndGetDistance();
     118             : 
     119      836788 : #pragma omp atomic
     120      836788 :             Aty[dataIndexForCurrentVoxel] += detectorValue * weight;
     121             : 
     122      836788 :             dataIndexForCurrentVoxel = traverse.getCurrentIndex();
     123      836788 :         }
     124       19001 :     }
     125             : 
     126             :     // ------------------------------------------
     127             :     // explicit template instantiation
     128             :     template class SiddonsMethod<float>;
     129             :     template class SiddonsMethod<double>;
     130             : } // namespace elsa