Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "LinearOperator.h"
       4             : #include "Geometry.h"
       5             : #include "BoundingBox.h"
       6             : #include "VolumeDescriptor.h"
       7             : #include "DetectorDescriptor.h"
       8             : #include "XrayProjector.h"
       9             : 
      10             : #include <vector>
      11             : #include <utility>
      12             : 
      13             : #include <Eigen/Geometry>
      14             : 
      15             : namespace elsa
      16             : {
      17             :     template <typename data_t = real_t>
      18             :     class JosephsMethod;
      19             : 
      20             :     template <typename data_t>
      21             :     struct XrayProjectorInnerTypes<JosephsMethod<data_t>> {
      22             :         using value_type = data_t;
      23             :         using forward_tag = any_projection_tag;
      24             :         using backward_tag = any_projection_tag;
      25             :     };
      26             : 
      27             :     /**
      28             :      * @brief Operator representing the discretized X-ray transform in 2d/3d using Joseph's method.
      29             :      *
      30             :      * @author Christoph Hahn - initial implementation
      31             :      * @author Maximilian Hornung - modularization
      32             :      * @author Nikola Dinev - fixes
      33             :      *
      34             :      * @tparam data_t data type for the domain and range of the operator, defaulting to real_t
      35             :      *
      36             :      * The volume is traversed along the rays as specified by the Geometry. For interior voxels
      37             :      * the sampling point is located in the middle of the two planes orthogonal to the main
      38             :      * direction of the ray. For boundary voxels the sampling point is located at the center of the
      39             :      * ray intersection with the voxel.
      40             :      *
      41             :      * The geometry is represented as a list of projection matrices (see class Geometry), one for
      42             :      * each acquisition pose.
      43             :      *
      44             :      * Two modes of interpolation are available:
      45             :      * NN (NearestNeighbours) takes the value of the pixel/voxel containing the point
      46             :      * LINEAR performs linear interpolation for the nearest 2 pixels (in 2D)
      47             :      * or the nearest 4 voxels (in 3D).
      48             :      *
      49             :      * Forward projection is accomplished using apply(), backward projection using applyAdjoint().
      50             :      * This projector is matched.
      51             :      */
      52             :     template <typename data_t>
      53             :     class JosephsMethod : public XrayProjector<JosephsMethod<data_t>>
      54             :     {
      55             :     public:
      56             :         /// Available interpolation modes
      57             :         enum class Interpolation { NN, LINEAR };
      58             : 
      59             :         using self_type = JosephsMethod<data_t>;
      60             :         using base_type = XrayProjector<self_type>;
      61             :         using value_type = typename base_type::value_type;
      62             :         using forward_tag = typename base_type::forward_tag;
      63             :         using backward_tag = typename base_type::backward_tag;
      64             : 
      65             :         /**
      66             :          * @brief Constructor for Joseph's traversal method.
      67             :          *
      68             :          * @param[in] domainDescriptor describing the domain of the operator (the volume)
      69             :          * @param[in] rangeDescriptor describing the range of the operator (the sinogram)
      70             :          * @param[in] interpolation enum specifying the interpolation mode
      71             :          *
      72             :          * The domain is expected to be 2 or 3 dimensional (volSizeX, volSizeY, [volSizeZ]),
      73             :          * the range is expected to be matching the domain (detSizeX, [detSizeY], acqPoses).
      74             :          */
      75             :         JosephsMethod(const VolumeDescriptor& domainDescriptor,
      76             :                       const DetectorDescriptor& rangeDescriptor,
      77             :                       Interpolation interpolation = Interpolation::LINEAR);
      78             : 
      79             :         /// default destructor
      80          59 :         ~JosephsMethod() = default;
      81             : 
      82             :     protected:
      83             :         /// default copy constructor, hidden from non-derived classes to prevent potential slicing
      84             :         JosephsMethod(const JosephsMethod<data_t>&) = default;
      85             : 
      86             :     private:
      87             :         /// apply Joseph's method (i.e. forward projection)
      88             :         void forward(const BoundingBox& aabb, const DataContainer<data_t>& x,
      89             :                      DataContainer<data_t>& Ax) const;
      90             : 
      91             :         /// apply the adjoint of Joseph's method (i.e. backward projection)
      92             :         void backward(const BoundingBox& aabb, const DataContainer<data_t>& y,
      93             :                       DataContainer<data_t>& Aty) const;
      94             : 
      95             :         /// implement the polymorphic clone operation
      96             :         JosephsMethod<data_t>* _cloneImpl() const;
      97             : 
      98             :         /// implement the polymorphic comparison operation
      99             :         bool _isEqual(const LinearOperator<data_t>& other) const;
     100             : 
     101             :         /// the interpolation mode
     102             :         Interpolation _interpolation;
     103             : 
     104             :         /// the traversal routine (for both apply/applyAdjoint)
     105             :         template <bool adjoint>
     106             :         void traverseVolume(const BoundingBox& aabb, const DataContainer<data_t>& vector,
     107             :                             DataContainer<data_t>& result) const;
     108             : 
     109             :         /**
     110             :          * @brief  Linear interpolation, works in any dimension
     111             :          *
     112             :          * @tparam adjoint true for backward projection, false for forward
     113             :          *
     114             :          * @param[in] vector the input DataContainer
     115             :          * @param[out] result DataContainer for results
     116             :          * @param[in] fractionals the fractional numbers used in the interpolation
     117             :          * @param[in] domainDim number of dimensions
     118             :          * @param[in] currentVoxel coordinates of voxel for interpolation
     119             :          * @param[in] intersection weighting for the interpolated values depending on the incidence
     120             :          * angle
     121             :          * @param[in] from index of the current vector position
     122             :          * @param[in] to index of the current result position
     123             :          * @param[in] mainDirection specifies the main direction of the ray
     124             :          */
     125             :         template <bool adjoint>
     126             :         void linear(const BoundingBox& aabb, const DataContainer<data_t>& vector,
     127             :                     DataContainer<data_t>& result, const RealVector_t& fractionals,
     128             :                     index_t domainDim, const IndexVector_t& currentVoxel, float intersection,
     129             :                     index_t from, index_t to, int mainDirection) const;
     130             : 
     131             :         /// lift from base class
     132             :         // using LinearOperator<data_t>::_domainDescriptor;
     133             :         // using LinearOperator<data_t>::_rangeDescriptor;
     134             : 
     135             :         friend class XrayProjector<self_type>;
     136             :     };
     137             : } // namespace elsa