LCOV - code coverage report
Current view: top level - projectors - JosephsMethod.cpp (source / functions) Hit Total Coverage
Test: test_coverage.info.cleaned Lines: 67 72 93.1 %
Date: 2022-05-27 02:48:28 Functions: 11 22 50.0 %

          Line data    Source code
       1             : #include "JosephsMethod.h"
       2             : #include "Timer.h"
       3             : #include "TraverseAABBJosephsMethod.h"
       4             : #include "Error.h"
       5             : #include "TypeCasts.hpp"
       6             : 
       7             : #include <type_traits>
       8             : 
       9             : namespace elsa
      10             : {
      11             :     template <typename data_t>
      12          63 :     JosephsMethod<data_t>::JosephsMethod(const VolumeDescriptor& domainDescriptor,
      13             :                                          const DetectorDescriptor& rangeDescriptor,
      14             :                                          Interpolation interpolation)
      15             :         : LinearOperator<data_t>(domainDescriptor, rangeDescriptor),
      16             :           _boundingBox{domainDescriptor.getNumberOfCoefficientsPerDimension()},
      17          63 :           _detectorDescriptor(static_cast<DetectorDescriptor&>(*_rangeDescriptor)),
      18          63 :           _volumeDescriptor(static_cast<VolumeDescriptor&>(*_domainDescriptor)),
      19          63 :           _interpolation{interpolation}
      20             :     {
      21          63 :         auto dim = _domainDescriptor->getNumberOfDimensions();
      22          63 :         if (dim != 2 && dim != 3) {
      23           0 :             throw InvalidArgumentError("JosephsMethod:only supporting 2d/3d operations");
      24             :         }
      25             : 
      26          63 :         if (dim != _rangeDescriptor->getNumberOfDimensions()) {
      27           0 :             throw InvalidArgumentError("JosephsMethod: domain and range dimension need to match");
      28             :         }
      29             : 
      30          63 :         if (_detectorDescriptor.getNumberOfGeometryPoses() == 0) {
      31           0 :             throw InvalidArgumentError("JosephsMethod: geometry list was empty");
      32             :         }
      33          63 :     }
      34             : 
      35             :     template <typename data_t>
      36          86 :     void JosephsMethod<data_t>::applyImpl(const DataContainer<data_t>& x,
      37             :                                           DataContainer<data_t>& Ax) const
      38             :     {
      39         258 :         Timer timeguard("JosephsMethod", "apply");
      40          86 :         traverseVolume<false>(x, Ax);
      41          86 :     }
      42             : 
      43             :     template <typename data_t>
      44          52 :     void JosephsMethod<data_t>::applyAdjointImpl(const DataContainer<data_t>& y,
      45             :                                                  DataContainer<data_t>& Aty) const
      46             :     {
      47         156 :         Timer timeguard("JosephsMethod", "applyAdjoint");
      48          52 :         traverseVolume<true>(y, Aty);
      49          52 :     }
      50             : 
      51             :     template <typename data_t>
      52          21 :     JosephsMethod<data_t>* JosephsMethod<data_t>::cloneImpl() const
      53             :     {
      54          21 :         return new JosephsMethod(_volumeDescriptor, _detectorDescriptor, _interpolation);
      55             :     }
      56             : 
      57             :     template <typename data_t>
      58           2 :     bool JosephsMethod<data_t>::isEqual(const LinearOperator<data_t>& other) const
      59             :     {
      60           2 :         if (!LinearOperator<data_t>::isEqual(other))
      61           0 :             return false;
      62             : 
      63           2 :         auto otherJM = downcast_safe<JosephsMethod>(&other);
      64           2 :         if (!otherJM)
      65           0 :             return false;
      66             : 
      67           2 :         return true;
      68             :     }
      69             : 
      70             :     template <typename data_t>
      71             :     template <bool adjoint>
      72         138 :     void JosephsMethod<data_t>::traverseVolume(const DataContainer<data_t>& vector,
      73             :                                                DataContainer<data_t>& result) const
      74             :     {
      75             :         if constexpr (adjoint)
      76          52 :             result = 0;
      77             : 
      78         138 :         const auto sizeOfRange = _rangeDescriptor->getNumberOfCoefficients();
      79         138 :         const auto rangeDim = _rangeDescriptor->getNumberOfDimensions();
      80             : 
      81             :         // iterate over all rays
      82         138 : #pragma omp parallel for
      83             :         for (index_t ir = 0; ir < sizeOfRange; ir++) {
      84             :             const auto ray = _detectorDescriptor.computeRayFromDetectorCoord(ir);
      85             : 
      86             :             // --> setup traversal algorithm
      87             :             TraverseAABBJosephsMethod traverse(_boundingBox, ray);
      88             : 
      89             :             if constexpr (!adjoint)
      90             :                 result[ir] = 0;
      91             : 
      92             :             // Make steps through the volume
      93             :             while (traverse.isInBoundingBox()) {
      94             : 
      95             :                 IndexVector_t currentVoxel = traverse.getCurrentVoxel();
      96             :                 float intersection = traverse.getIntersectionLength();
      97             : 
      98             :                 // to avoid code duplicates for apply and applyAdjoint
      99      176080 :                 auto [to, from] = [&]() {
     100      528240 :                     const auto curVoxIndex =
     101      176080 :                         _domainDescriptor->getIndexFromCoordinate(currentVoxel);
     102      176080 :                     return adjoint ? std::pair{curVoxIndex, ir} : std::pair{ir, curVoxIndex};
     103      176080 :                 }();
     104             : 
     105             :                 // #dfrank: This supresses a clang-tidy warning:
     106             :                 // "error: reference to local binding 'to' declared in enclosing
     107             :                 // context [clang-diagnostic-error]", which seems to be based on
     108             :                 // clang-8, and a bug in the standard, as structured bindings are not
     109             :                 // "real" variables, but only references. I'm not sure why a warning is
     110             :                 // generated though
     111             :                 auto tmpTo = to;
     112             :                 auto tmpFrom = from;
     113             : 
     114             :                 switch (_interpolation) {
     115             :                     case Interpolation::LINEAR:
     116             :                         linear<adjoint>(vector, result, traverse.getFractionals(), rangeDim,
     117             :                                         currentVoxel, intersection, from, to,
     118             :                                         traverse.getIgnoreDirection());
     119             :                         break;
     120             :                     case Interpolation::NN:
     121             :                         if constexpr (adjoint) {
     122             : #pragma omp atomic
     123             :                             // NOLINTNEXTLINE
     124             :                             result[tmpTo] += intersection * vector[tmpFrom];
     125             :                         } else {
     126             :                             result[tmpTo] += intersection * vector[tmpFrom];
     127             :                         }
     128             :                         break;
     129             :                 }
     130             : 
     131             :                 // update Traverse
     132             :                 traverse.updateTraverse();
     133             :             }
     134             :         }
     135         138 :     }
     136             : 
     137             :     template <typename data_t>
     138             :     template <bool adjoint>
     139      176080 :     void JosephsMethod<data_t>::linear(const DataContainer<data_t>& vector,
     140             :                                        DataContainer<data_t>& result,
     141             :                                        const RealVector_t& fractionals, index_t domainDim,
     142             :                                        const IndexVector_t& currentVoxel, float intersection,
     143             :                                        index_t from, index_t to, int mainDirection) const
     144             :     {
     145      176080 :         data_t weight = intersection;
     146      352160 :         IndexVector_t interpol = currentVoxel;
     147             : 
     148             :         // handle diagonal if 3D
     149      176080 :         if (domainDim == 3) {
     150         456 :             for (int i = 0; i < domainDim; i++) {
     151         342 :                 if (i != mainDirection) {
     152         228 :                     weight *= std::abs(fractionals(i));
     153         228 :                     interpol(i) += (fractionals(i) < 0.0) ? -1 : 1;
     154             :                     // mirror values at border if outside the volume
     155         228 :                     auto interpolVal = static_cast<real_t>(interpol(i));
     156         228 :                     if (interpolVal < _boundingBox._min(i) || interpolVal > _boundingBox._max(i))
     157           7 :                         interpol(i) = static_cast<index_t>(_boundingBox._min(i));
     158         221 :                     else if (interpolVal == _boundingBox._max(i))
     159          13 :                         interpol(i) = static_cast<index_t>(_boundingBox._max(i)) - 1;
     160             :                 }
     161             :             }
     162             :             if constexpr (adjoint) {
     163             : #pragma omp atomic
     164          43 :                 result(interpol) += weight * vector[from];
     165             :             } else {
     166          71 :                 result[to] += weight * vector(interpol);
     167             :             }
     168             :         }
     169             : 
     170             :         // handle current voxel
     171      176080 :         weight = intersection * (1 - fractionals.array().abs()).prod()
     172      176080 :                  / (1 - std::abs(fractionals(mainDirection)));
     173             :         if constexpr (adjoint) {
     174             : #pragma omp atomic
     175       87983 :             result[to] += weight * vector[from];
     176             :         } else {
     177       88105 :             result[to] += weight * vector[from];
     178             :         }
     179             : 
     180             :         // handle neighbors not along the main direction
     181      528354 :         for (int i = 0; i < domainDim; i++) {
     182      352274 :             if (i != mainDirection) {
     183      176194 :                 data_t weightn = weight * std::abs(fractionals(i)) / (1 - std::abs(fractionals(i)));
     184      176194 :                 interpol = currentVoxel;
     185      176194 :                 interpol(i) += (fractionals(i) < 0.0) ? -1 : 1;
     186             : 
     187             :                 // mirror values at border if outside the volume
     188      176194 :                 auto interpolVal = static_cast<real_t>(interpol(i));
     189      176194 :                 if (interpolVal < _boundingBox._min(i) || interpolVal > _boundingBox._max(i))
     190        2951 :                     interpol(i) = static_cast<index_t>(_boundingBox._min(i));
     191      173243 :                 else if (interpolVal == _boundingBox._max(i))
     192        2917 :                     interpol(i) = static_cast<index_t>(_boundingBox._max(i)) - 1;
     193             : 
     194             :                 if constexpr (adjoint) {
     195             : #pragma omp atomic
     196       88045 :                     result(interpol) += weightn * vector[from];
     197             :                 } else {
     198       88176 :                     result[to] += weightn * vector(interpol);
     199             :                 }
     200             :             }
     201             :         }
     202      176080 :     }
     203             : 
     204             :     // ------------------------------------------
     205             :     // explicit template instantiation
     206             :     template class JosephsMethod<float>;
     207             :     template class JosephsMethod<double>;
     208             : 
     209             : } // namespace elsa

Generated by: LCOV version 1.15