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