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 0 : SubsetSampler<DetectorDescriptor_t, data_t>::SubsetSampler(
11 : const VolumeDescriptor& volumeDescriptor, const DetectorDescriptor_t& detectorDescriptor,
12 : index_t nSubsets, SamplingStrategy samplingStrategy)
13 0 : : _indexMapping(static_cast<std::size_t>(nSubsets)),
14 0 : _volumeDescriptor(volumeDescriptor),
15 0 : _fullDetectorDescriptor(detectorDescriptor),
16 0 : _nSubsets{nSubsets}
17 : {
18 0 : 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 0 : const auto numCoeffsPerDim = detectorDescriptor.getNumberOfCoefficientsPerDimension();
25 0 : const index_t nDimensions = detectorDescriptor.getNumberOfDimensions();
26 0 : const auto numElements = numCoeffsPerDim[nDimensions - 1];
27 0 : if (samplingStrategy == SamplingStrategy::ROUND_ROBIN) {
28 0 : std::vector<index_t> indices(static_cast<std::size_t>(numElements));
29 0 : std::iota(indices.begin(), indices.end(), 0);
30 0 : _indexMapping = splitRoundRobin(indices, _nSubsets);
31 0 : } else if (samplingStrategy == SamplingStrategy::ROTATIONAL_CLUSTERING) {
32 0 : _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 0 : for (const auto& blockIndices : _indexMapping) {
39 0 : std::vector<Geometry> geometry;
40 0 : for (auto index : blockIndices) {
41 0 : auto geo = detectorDescriptor.getGeometryAt(index);
42 0 : if (geo.has_value()) {
43 0 : geometry.emplace_back(*geo);
44 : }
45 : }
46 0 : IndexVector_t numOfCoeffsPerDim =
47 : detectorDescriptor.getNumberOfCoefficientsPerDimension();
48 0 : numOfCoeffsPerDim[numOfCoeffsPerDim.size() - 1] =
49 0 : static_cast<index_t>(blockIndices.size());
50 :
51 0 : _detectorDescriptors.emplace_back(DetectorDescriptor_t(numOfCoeffsPerDim, geometry));
52 : }
53 0 : }
54 :
55 : template <typename DetectorDescriptor_t, typename data_t>
56 0 : std::vector<std::vector<index_t>> SubsetSampler<DetectorDescriptor_t, data_t>::splitRoundRobin(
57 : const std::vector<index_t>& indices, index_t nSubsets)
58 : {
59 0 : std::vector<std::vector<index_t>> subsetIndices(static_cast<std::size_t>(nSubsets));
60 :
61 : // determine the mapping of indices to subsets
62 0 : for (std::size_t i = 0; i < indices.size(); ++i) {
63 0 : const auto subset = i % static_cast<std::size_t>(nSubsets);
64 0 : subsetIndices[subset].template emplace_back(indices[i]);
65 : }
66 :
67 0 : return subsetIndices;
68 0 : }
69 :
70 : template <typename DetectorDescriptor_t, typename data_t>
71 : std::vector<std::vector<index_t>>
72 0 : SubsetSampler<DetectorDescriptor_t, data_t>::splitRotationalClustering(
73 : const DetectorDescriptor_t& detectorDescriptor, index_t nSubsets)
74 : {
75 :
76 0 : const auto numCoeffsPerDim = detectorDescriptor.getNumberOfCoefficientsPerDimension();
77 0 : const index_t nDimensions = detectorDescriptor.getNumberOfDimensions();
78 0 : const auto numElements = numCoeffsPerDim[nDimensions - 1];
79 0 : std::vector<index_t> indices(static_cast<std::size_t>(numElements));
80 0 : std::iota(indices.begin(), indices.end(), 0);
81 0 : const auto geometry = detectorDescriptor.getGeometry();
82 :
83 : // angle between two rotation matrices used as a distance measure
84 0 : auto dist = [nDimensions](auto& g1, auto& g2) {
85 0 : const auto& r1 = g1.getRotationMatrix();
86 0 : const auto& r2 = g2.getRotationMatrix();
87 0 : auto product = r1 * r2.transpose();
88 0 : if (nDimensions == 2) { // the 2D case
89 0 : return static_cast<double>(std::atan2(product(1, 0), product(0, 0)));
90 : } else { // the 3D case
91 0 : return std::acos((product.trace() - 1.0) / 2.0);
92 : }
93 : };
94 :
95 0 : const auto first = geometry.front();
96 0 : std::sort(std::begin(indices), std::end(indices),
97 0 : [dist, first, &geometry](auto lhs, auto rhs) {
98 0 : return dist(first, geometry[static_cast<std::size_t>(lhs)])
99 0 : < dist(first, geometry[static_cast<std::size_t>(rhs)]);
100 : });
101 :
102 0 : return splitRoundRobin(indices, nSubsets);
103 0 : }
104 :
105 : template <typename DetectorDescriptor_t, typename data_t>
106 0 : SubsetSampler<DetectorDescriptor_t, data_t>::SubsetSampler(
107 : const SubsetSampler<DetectorDescriptor_t, data_t>& other)
108 0 : : _indexMapping{other._indexMapping},
109 0 : _volumeDescriptor(other._volumeDescriptor),
110 0 : _fullDetectorDescriptor(other._fullDetectorDescriptor),
111 0 : _nSubsets{other._nSubsets}
112 : {
113 0 : for (const auto& detectorDescriptor : other._detectorDescriptors) {
114 0 : _detectorDescriptors.emplace_back(detectorDescriptor);
115 : }
116 0 : }
117 :
118 : template <typename DetectorDescriptor_t, typename data_t>
119 0 : 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 0 : IndexVector_t slicesInBlock(_indexMapping.size());
124 0 : for (unsigned long i = 0; i < _indexMapping.size(); i++) {
125 0 : slicesInBlock[static_cast<index_t>(i)] = static_cast<index_t>(_indexMapping[i].size());
126 : }
127 0 : PartitionDescriptor dataDescriptor(sinogram.getDataDescriptor(), slicesInBlock);
128 :
129 0 : const auto numCoeffsPerDim =
130 0 : sinogram.getDataDescriptor().getNumberOfCoefficientsPerDimension();
131 0 : auto partitionedData = DataContainer<data_t>(dataDescriptor, sinogram.getDataHandlerType());
132 :
133 : // resort the actual measurement partitionedData
134 0 : index_t coeffsPerRow = numCoeffsPerDim.head(numCoeffsPerDim.size() - 1).prod();
135 0 : for (index_t i = 0; i < _nSubsets; i++) {
136 : // the indices of the partitionedData rows belonging to this subset
137 0 : std::vector<index_t> indices = _indexMapping[static_cast<std::size_t>(i)];
138 :
139 0 : auto block = partitionedData.getBlock(i);
140 :
141 0 : for (std::size_t j = 0; j < indices.size(); j++) {
142 0 : for (int k = 0; k < coeffsPerRow; k++) {
143 0 : block[static_cast<index_t>(j) * coeffsPerRow + k] =
144 0 : sinogram[indices[j] * coeffsPerRow + k];
145 : }
146 : }
147 : }
148 0 : return partitionedData;
149 0 : }
150 :
151 : template <typename DetectorDescriptor_t, typename data_t>
152 : SubsetSampler<DetectorDescriptor_t, data_t>*
153 0 : SubsetSampler<DetectorDescriptor_t, data_t>::cloneImpl() const
154 : {
155 0 : return new SubsetSampler<DetectorDescriptor_t, data_t>(*this);
156 : }
157 :
158 : template <typename DetectorDescriptor_t, typename data_t>
159 0 : bool SubsetSampler<DetectorDescriptor_t, data_t>::isEqual(
160 : const SubsetSampler<DetectorDescriptor_t, data_t>& other) const
161 : {
162 0 : if (typeid(*this) != typeid(other))
163 0 : return false;
164 :
165 0 : if (_indexMapping != other._indexMapping)
166 0 : return false;
167 :
168 0 : if (_volumeDescriptor != other._volumeDescriptor)
169 0 : return false;
170 :
171 0 : if (_fullDetectorDescriptor != other._fullDetectorDescriptor)
172 0 : return false;
173 :
174 0 : 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 0 : 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
|
