Line data    Source code

       1             : #include <set>
       2             : #include <numeric>
       3             : #include "SubsetSampler.h"
       4             : #include "PartitionDescriptor.h"
       5             : #include "SiddonsMethod.h"
       6             : 
       7             : namespace elsa
       8             : {
       9             :     template <typename DetectorDescriptor_t, typename data_t>
      10             :     SubsetSampler<DetectorDescriptor_t, data_t>::SubsetSampler(
      11             :         const VolumeDescriptor& volumeDescriptor, const DetectorDescriptor_t& detectorDescriptor,
      12             :         index_t nSubsets, SamplingStrategy samplingStrategy)
      13             :         : _indexMapping(static_cast<std::size_t>(nSubsets)),
      14             :           _volumeDescriptor(volumeDescriptor),
      15             :           _fullDetectorDescriptor(detectorDescriptor),
      16             :           _nSubsets{nSubsets}
      17          14 :     {
      18          14 :         if (nSubsets <= 1) {
      19           0 :             throw std::invalid_argument("SubsetSampler: nSubsets must be >= 2");
      20           0 :         }
      21             : 
      22             :         // the mapping of data indices to subsets
      23             : 
      24          14 :         const auto numCoeffsPerDim = detectorDescriptor.getNumberOfCoefficientsPerDimension();
      25          14 :         const index_t nDimensions = detectorDescriptor.getNumberOfDimensions();
      26          14 :         const auto numElements = numCoeffsPerDim[nDimensions - 1];
      27          14 :         if (samplingStrategy == SamplingStrategy::ROUND_ROBIN) {
      28           7 :             std::vector<index_t> indices(static_cast<std::size_t>(numElements));
      29           7 :             std::iota(indices.begin(), indices.end(), 0);
      30           7 :             _indexMapping = splitRoundRobin(indices, _nSubsets);
      31           7 :         } else if (samplingStrategy == SamplingStrategy::ROTATIONAL_CLUSTERING) {
      32           7 :             _indexMapping = splitRotationalClustering(detectorDescriptor, _nSubsets);
      33           7 :         } else {
      34           0 :             throw std::invalid_argument("SubsetSampler: unsupported sampling strategy");
      35           0 :         }
      36             : 
      37             :         // create the detector descriptors that correspond to each subset
      38          56 :         for (const auto& blockIndices : _indexMapping) {
      39          56 :             std::vector<Geometry> geometry;
      40         256 :             for (auto index : blockIndices) {
      41         256 :                 auto geo = detectorDescriptor.getGeometryAt(index);
      42         256 :                 if (geo.has_value()) {
      43         256 :                     geometry.emplace_back(*geo);
      44         256 :                 }
      45         256 :             }
      46          56 :             IndexVector_t numOfCoeffsPerDim =
      47          56 :                 detectorDescriptor.getNumberOfCoefficientsPerDimension();
      48          56 :             numOfCoeffsPerDim[numOfCoeffsPerDim.size() - 1] =
      49          56 :                 static_cast<index_t>(blockIndices.size());
      50             : 
      51          56 :             _detectorDescriptors.emplace_back(DetectorDescriptor_t(numOfCoeffsPerDim, geometry));
      52          56 :         }
      53          14 :     }
      54             : 
      55             :     template <typename DetectorDescriptor_t, typename data_t>
      56             :     std::vector<std::vector<index_t>> SubsetSampler<DetectorDescriptor_t, data_t>::splitRoundRobin(
      57             :         const std::vector<index_t>& indices, index_t nSubsets)
      58          20 :     {
      59          20 :         std::vector<std::vector<index_t>> subsetIndices(static_cast<std::size_t>(nSubsets));
      60             : 
      61             :         // determine the mapping of indices to subsets
      62         468 :         for (std::size_t i = 0; i < indices.size(); ++i) {
      63         448 :             const auto subset = i % static_cast<std::size_t>(nSubsets);
      64         448 :             subsetIndices[subset].template emplace_back(indices[i]);
      65         448 :         }
      66             : 
      67          20 :         return subsetIndices;
      68          20 :     }
      69             : 
      70             :     template <typename DetectorDescriptor_t, typename data_t>
      71             :     std::vector<std::vector<index_t>>
      72             :         SubsetSampler<DetectorDescriptor_t, data_t>::splitRotationalClustering(
      73             :             const DetectorDescriptor_t& detectorDescriptor, index_t nSubsets)
      74           9 :     {
      75             : 
      76           9 :         const auto numCoeffsPerDim = detectorDescriptor.getNumberOfCoefficientsPerDimension();
      77           9 :         const index_t nDimensions = detectorDescriptor.getNumberOfDimensions();
      78           9 :         const auto numElements = numCoeffsPerDim[nDimensions - 1];
      79           9 :         std::vector<index_t> indices(static_cast<std::size_t>(numElements));
      80           9 :         std::iota(indices.begin(), indices.end(), 0);
      81           9 :         const auto geometry = detectorDescriptor.getGeometry();
      82             : 
      83             :         // angle between two rotation matrices used as a distance measure
      84        2026 :         auto dist = [nDimensions](auto& g1, auto& g2) {
      85        2026 :             const auto& r1 = g1.getRotationMatrix();
      86        2026 :             const auto& r2 = g2.getRotationMatrix();
      87        2026 :             auto product = r1 * r2.transpose();
      88        2026 :             if (nDimensions == 2) { // the 2D case
      89        1304 :                 return static_cast<double>(std::atan2(product(1, 0), product(0, 0)));
      90        1304 :             } else { // the 3D case
      91         722 :                 return std::acos((product.trace() - 1.0) / 2.0);
      92         722 :             }
      93        2026 :         };
      94             : 
      95           9 :         const auto first = geometry.front();
      96           9 :         std::sort(std::begin(indices), std::end(indices),
      97        1013 :                   [dist, first, &geometry](auto lhs, auto rhs) {
      98        1013 :                       return dist(first, geometry[static_cast<std::size_t>(lhs)])
      99        1013 :                              < dist(first, geometry[static_cast<std::size_t>(rhs)]);
     100        1013 :                   });
     101             : 
     102           9 :         return splitRoundRobin(indices, nSubsets);
     103           9 :     }
     104             : 
     105             :     template <typename DetectorDescriptor_t, typename data_t>
     106             :     SubsetSampler<DetectorDescriptor_t, data_t>::SubsetSampler(
     107             :         const SubsetSampler<DetectorDescriptor_t, data_t>& other)
     108             :         : _indexMapping{other._indexMapping},
     109             :           _volumeDescriptor(other._volumeDescriptor),
     110             :           _fullDetectorDescriptor(other._fullDetectorDescriptor),
     111             :           _nSubsets{other._nSubsets}
     112           4 :     {
     113          16 :         for (const auto& detectorDescriptor : other._detectorDescriptors) {
     114          16 :             _detectorDescriptors.emplace_back(detectorDescriptor);
     115          16 :         }
     116           4 :     }
     117             : 
     118             :     template <typename DetectorDescriptor_t, typename data_t>
     119             :     DataContainer<data_t> SubsetSampler<DetectorDescriptor_t, data_t>::getPartitionedData(
     120             :         const DataContainer<data_t>& sinogram)
     121          10 :     {
     122             :         // save the number of entries per subset
     123          10 :         IndexVector_t slicesInBlock(_indexMapping.size());
     124          50 :         for (unsigned long i = 0; i < _indexMapping.size(); i++) {
     125          40 :             slicesInBlock[static_cast<index_t>(i)] = static_cast<index_t>(_indexMapping[i].size());
     126          40 :         }
     127          10 :         PartitionDescriptor dataDescriptor(sinogram.getDataDescriptor(), slicesInBlock);
     128             : 
     129          10 :         const auto numCoeffsPerDim =
     130          10 :             sinogram.getDataDescriptor().getNumberOfCoefficientsPerDimension();
     131          10 :         auto partitionedData = DataContainer<data_t>(dataDescriptor, sinogram.getDataHandlerType());
     132             : 
     133             :         // resort the actual measurement partitionedData
     134          10 :         index_t coeffsPerRow = numCoeffsPerDim.head(numCoeffsPerDim.size() - 1).prod();
     135          50 :         for (index_t i = 0; i < _nSubsets; i++) {
     136             :             // the indices of the partitionedData rows belonging to this subset
     137          40 :             std::vector<index_t> indices = _indexMapping[static_cast<std::size_t>(i)];
     138             : 
     139          40 :             auto block = partitionedData.getBlock(i);
     140             : 
     141         224 :             for (std::size_t j = 0; j < indices.size(); j++) {
     142       10568 :                 for (int k = 0; k < coeffsPerRow; k++) {
     143       10384 :                     block[static_cast<index_t>(j) * coeffsPerRow + k] =
     144       10384 :                         sinogram[indices[j] * coeffsPerRow + k];
     145       10384 :                 }
     146         184 :             }
     147          40 :         }
     148          10 :         return partitionedData;
     149          10 :     }
     150             : 
     151             :     template <typename DetectorDescriptor_t, typename data_t>
     152             :     SubsetSampler<DetectorDescriptor_t, data_t>*
     153             :         SubsetSampler<DetectorDescriptor_t, data_t>::cloneImpl() const
     154           4 :     {
     155           4 :         return new SubsetSampler<DetectorDescriptor_t, data_t>(*this);
     156           4 :     }
     157             : 
     158             :     template <typename DetectorDescriptor_t, typename data_t>
     159             :     bool SubsetSampler<DetectorDescriptor_t, data_t>::isEqual(
     160             :         const SubsetSampler<DetectorDescriptor_t, data_t>& other) const
     161           4 :     {
     162           4 :         if (typeid(*this) != typeid(other))
     163           0 :             return false;
     164             : 
     165           4 :         if (_indexMapping != other._indexMapping)
     166           0 :             return false;
     167             : 
     168           4 :         if (_volumeDescriptor != other._volumeDescriptor)
     169           0 :             return false;
     170             : 
     171           4 :         if (_fullDetectorDescriptor != other._fullDetectorDescriptor)
     172           0 :             return false;
     173             : 
     174           4 :         if (_nSubsets != other._nSubsets)
     175           0 :             return false;
     176             : 
     177             :         // we do not need to check if the vector of detector descriptors is equal as this is
     178             :         // implied by the equality of _data in combination with the full detector descriptor
     179             : 
     180           4 :         return true;
     181           4 :     }
     182             : 
     183             :     // ------------------------------------------
     184             :     // explicit template instantiation
     185             :     template class SubsetSampler<PlanarDetectorDescriptor, float>;
     186             :     template class SubsetSampler<PlanarDetectorDescriptor, double>;
     187             :     template class SubsetSampler<PlanarDetectorDescriptor, complex<float>>;
     188             :     template class SubsetSampler<PlanarDetectorDescriptor, complex<double>>;
     189             : } // namespace elsa