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