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 14 : 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 14 : _nSubsets{nSubsets}
17 : {
18 14 : if (nSubsets <= 1) {
19 0 : throw std::invalid_argument("SubsetSampler: nSubsets must be >= 2");
20 : }
21 :
22 : // the mapping of data indices to subsets
23 :
24 28 : 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 : } else {
34 0 : throw std::invalid_argument("SubsetSampler: unsupported sampling strategy");
35 : }
36 :
37 : // create the detector descriptors that correspond to each subset
38 70 : for (const auto& blockIndices : _indexMapping) {
39 112 : std::vector<Geometry> geometry;
40 312 : for (auto index : blockIndices) {
41 512 : auto geo = detectorDescriptor.getGeometryAt(index);
42 256 : if (geo.has_value()) {
43 256 : geometry.emplace_back(*geo);
44 : }
45 : }
46 56 : IndexVector_t numOfCoeffsPerDim =
47 : 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 : }
53 14 : }
54 :
55 : template <typename DetectorDescriptor_t, typename data_t>
56 20 : std::vector<std::vector<index_t>> SubsetSampler<DetectorDescriptor_t, data_t>::splitRoundRobin(
57 : const std::vector<index_t>& indices, index_t nSubsets)
58 : {
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 : }
66 :
67 20 : return subsetIndices;
68 : }
69 :
70 : template <typename DetectorDescriptor_t, typename data_t>
71 : std::vector<std::vector<index_t>>
72 9 : SubsetSampler<DetectorDescriptor_t, data_t>::splitRotationalClustering(
73 : const DetectorDescriptor_t& detectorDescriptor, index_t nSubsets)
74 : {
75 :
76 18 : const auto numCoeffsPerDim = detectorDescriptor.getNumberOfCoefficientsPerDimension();
77 9 : const index_t nDimensions = detectorDescriptor.getNumberOfDimensions();
78 9 : const auto numElements = numCoeffsPerDim[nDimensions - 1];
79 18 : std::vector<index_t> indices(static_cast<std::size_t>(numElements));
80 9 : std::iota(indices.begin(), indices.end(), 0);
81 18 : const auto geometry = detectorDescriptor.getGeometry();
82 :
83 : // angle between two rotation matrices used as a distance measure
84 4061 : 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 : } else { // the 3D case
91 722 : return std::acos((product.trace() - 1.0) / 2.0);
92 : }
93 : };
94 :
95 18 : const auto first = geometry.front();
96 2035 : 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 : });
101 :
102 18 : return splitRoundRobin(indices, nSubsets);
103 : }
104 :
105 : template <typename DetectorDescriptor_t, typename data_t>
106 4 : SubsetSampler<DetectorDescriptor_t, data_t>::SubsetSampler(
107 : const SubsetSampler<DetectorDescriptor_t, data_t>& other)
108 4 : : _indexMapping{other._indexMapping},
109 4 : _volumeDescriptor(other._volumeDescriptor),
110 4 : _fullDetectorDescriptor(other._fullDetectorDescriptor),
111 4 : _nSubsets{other._nSubsets}
112 : {
113 20 : for (const auto& detectorDescriptor : other._detectorDescriptors) {
114 16 : _detectorDescriptors.emplace_back(detectorDescriptor);
115 : }
116 4 : }
117 :
118 : template <typename DetectorDescriptor_t, typename data_t>
119 10 : DataContainer<data_t> SubsetSampler<DetectorDescriptor_t, data_t>::getPartitionedData(
120 : const DataContainer<data_t>& sinogram)
121 : {
122 : // save the number of entries per subset
123 20 : 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 : }
127 20 : 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 80 : std::vector<index_t> indices = _indexMapping[static_cast<std::size_t>(i)];
138 :
139 80 : 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 : }
146 : }
147 : }
148 20 : return partitionedData;
149 : }
150 :
151 : template <typename DetectorDescriptor_t, typename data_t>
152 : SubsetSampler<DetectorDescriptor_t, data_t>*
153 4 : SubsetSampler<DetectorDescriptor_t, data_t>::cloneImpl() const
154 : {
155 4 : return new SubsetSampler<DetectorDescriptor_t, data_t>(*this);
156 : }
157 :
158 : template <typename DetectorDescriptor_t, typename data_t>
159 4 : bool SubsetSampler<DetectorDescriptor_t, data_t>::isEqual(
160 : const SubsetSampler<DetectorDescriptor_t, data_t>& other) const
161 : {
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 : }
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
|
