LCOV - code coverage report
Current view: top level - projectors - SiddonsMethod.cpp (source / functions) Hit Total Coverage
Test: test_coverage.info.cleaned Lines: 28 33 84.8 %
Date: 2022-05-27 02:48:28 Functions: 7 14 50.0 %

          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         347 :     SiddonsMethod<data_t>::SiddonsMethod(const VolumeDescriptor& domainDescriptor,
      13             :                                          const DetectorDescriptor& rangeDescriptor)
      14             :         : LinearOperator<data_t>(domainDescriptor, rangeDescriptor),
      15             :           _boundingBox(domainDescriptor.getNumberOfCoefficientsPerDimension()),
      16         347 :           _detectorDescriptor(static_cast<DetectorDescriptor&>(*_rangeDescriptor)),
      17         347 :           _volumeDescriptor(static_cast<VolumeDescriptor&>(*_domainDescriptor))
      18             :     {
      19         347 :         auto dim = _domainDescriptor->getNumberOfDimensions();
      20         347 :         if (dim != _rangeDescriptor->getNumberOfDimensions()) {
      21           0 :             throw LogicError("SiddonsMethod: domain and range dimension need to match");
      22             :         }
      23             : 
      24         347 :         if (dim != 2 && dim != 3) {
      25           0 :             throw LogicError("SiddonsMethod: only supporting 2d/3d operations");
      26             :         }
      27             : 
      28         347 :         if (_detectorDescriptor.getNumberOfGeometryPoses() == 0) {
      29           0 :             throw LogicError("SiddonsMethod: geometry list was empty");
      30             :         }
      31         347 :     }
      32             : 
      33             :     template <typename data_t>
      34         266 :     void SiddonsMethod<data_t>::applyImpl(const DataContainer<data_t>& x,
      35             :                                           DataContainer<data_t>& Ax) const
      36             :     {
      37         798 :         Timer t("SiddonsMethod", "apply");
      38         266 :         traverseVolume<false>(x, Ax);
      39         266 :     }
      40             : 
      41             :     template <typename data_t>
      42         218 :     void SiddonsMethod<data_t>::applyAdjointImpl(const DataContainer<data_t>& y,
      43             :                                                  DataContainer<data_t>& Aty) const
      44             :     {
      45         654 :         Timer t("SiddonsMethod", "applyAdjoint");
      46         218 :         traverseVolume<true>(y, Aty);
      47         218 :     }
      48             : 
      49             :     template <typename data_t>
      50         266 :     SiddonsMethod<data_t>* SiddonsMethod<data_t>::cloneImpl() const
      51             :     {
      52         266 :         return new SiddonsMethod(_volumeDescriptor, _detectorDescriptor);
      53             :     }
      54             : 
      55             :     template <typename data_t>
      56           5 :     bool SiddonsMethod<data_t>::isEqual(const LinearOperator<data_t>& other) const
      57             :     {
      58           5 :         if (!LinearOperator<data_t>::isEqual(other))
      59           0 :             return false;
      60             : 
      61           5 :         auto otherSM = downcast_safe<SiddonsMethod>(&other);
      62           5 :         if (!otherSM)
      63           0 :             return false;
      64             : 
      65           5 :         return true;
      66             :     }
      67             : 
      68             :     template <typename data_t>
      69             :     template <bool adjoint>
      70         484 :     void SiddonsMethod<data_t>::traverseVolume(const DataContainer<data_t>& vector,
      71             :                                                DataContainer<data_t>& result) const
      72             :     {
      73         484 :         const index_t maxIterations = adjoint ? vector.getSize() : result.getSize();
      74             : 
      75             :         if constexpr (adjoint) {
      76         218 :             result = 0; // initialize volume to 0, because we are not going to hit every voxel!
      77             :         }
      78             : 
      79             :         // --> loop either over every voxel that should  updated or every detector
      80             :         // cell that should be calculated
      81         484 : #pragma omp parallel for
      82             :         for (index_t rangeIndex = 0; rangeIndex < maxIterations; ++rangeIndex) {
      83             :             // --> get the current ray to the detector center
      84             :             const auto ray = _detectorDescriptor.computeRayFromDetectorCoord(rangeIndex);
      85             : 
      86             :             // --> setup traversal algorithm
      87             :             TraverseAABB traverse(_boundingBox, ray);
      88             : 
      89             :             if constexpr (!adjoint)
      90             :                 result[rangeIndex] = 0;
      91             : 
      92             :             // --> initial index to access the data vector
      93             :             auto dataIndexForCurrentVoxel =
      94             :                 _domainDescriptor->getIndexFromCoordinate(traverse.getCurrentVoxel());
      95             : 
      96             :             while (traverse.isInBoundingBox()) {
      97             : 
      98             :                 auto weight = traverse.updateTraverseAndGetDistance();
      99             :                 // --> update result depending on the operation performed
     100             :                 if constexpr (adjoint)
     101             : #pragma omp atomic
     102             :                     result[dataIndexForCurrentVoxel] += vector[rangeIndex] * weight;
     103             :                 else
     104             :                     result[rangeIndex] += vector[dataIndexForCurrentVoxel] * weight;
     105             : 
     106             :                 dataIndexForCurrentVoxel =
     107             :                     _domainDescriptor->getIndexFromCoordinate(traverse.getCurrentVoxel());
     108             :             }
     109             :         }
     110         484 :     }
     111             : 
     112             :     // ------------------------------------------
     113             :     // explicit template instantiation
     114             :     template class SiddonsMethod<float>;
     115             :     template class SiddonsMethod<double>;
     116             : } // namespace elsa

Generated by: LCOV version 1.15