Line data    Source code

       1             : #include "FBP.h"
       2             : #include "DataContainer.h"
       3             : #include "Filter.h"
       4             : #include "FourierTransform.h"
       5             : #include "TypeCasts.hpp"
       6             : #include <Eigen/src/Core/util/Constants.h>
       7             : 
       8             : namespace elsa
       9             : {
      10             :     template <typename data_t>
      11             :     FBP<data_t>::FBP(const LinearOperator<data_t>& P, const Filter<data_t>& g)
      12             :         : projector_{P.clone()}, filter_{downcast<Filter<data_t>>(g.clone())}
      13           0 :     {
      14           0 :     }
      15             : 
      16             :     template <typename data_t>
      17             :     FBP<data_t>::FBP(const LinearOperator<data_t>& P, std::unique_ptr<Filter<data_t>> g)
      18             :         : projector_{P.clone()}, filter_{std::move(g)}
      19           1 :     {
      20           1 :     }
      21             : 
      22             :     template <typename data_t>
      23             :     DataContainer<data_t> FBP<data_t>::apply(const DataContainer<data_t>& sinogram) const
      24           1 :     {
      25           1 :         auto& descriptor = sinogram.getDataDescriptor();
      26           1 :         auto dim = descriptor.getNumberOfDimensions();
      27             : 
      28           1 :         if (dim != 2) {
      29           0 :             throw InvalidArgumentError("FBP.apply:: only handles 2D data");
      30           0 :         }
      31             : 
      32           1 :         auto filtered = empty<data_t>(descriptor);
      33           1 :         auto numSlices = descriptor.getNumberOfCoefficientsPerDimension()[dim - 1];
      34             : 
      35           1 :         PartitionDescriptor sliceDescriptor{descriptor, numSlices};
      36           1 :         FourierTransform<complex<data_t>> F{sliceDescriptor.getDescriptorOfBlock(0)};
      37             : 
      38         201 :         for (auto i = 0; i < numSlices; i++) {
      39         200 :             auto spectrum = F.apply(sinogram.slice(i).asComplex());
      40         200 :             auto filteredRow = F.applyAdjoint(filter_->apply(spectrum));
      41         200 :             filtered.slice(i) = elsa::real(filteredRow);
      42         200 :         }
      43             : 
      44           1 :         auto backprojected = projector_->applyAdjoint(filtered);
      45           1 :         return backprojected * pi_t / 2 / numSlices; // This normalization is necessary because the
      46             :                                                      // projectors are not normalized
      47             :                                                      // TODO: Compare normalization with literature
      48           1 :     }
      49             : 
      50             :     // ------------------------------------------
      51             :     // explicit template instantiation
      52             :     template class FBP<float>;
      53             :     template class FBP<double>;
      54             : 
      55             : } // namespace elsa