Line data Source code
1 : #include "TypeCasts.hpp"
2 : #include "elsaDefines.h"
3 : #include "Error.h"
4 :
5 : #include <Eigen/Core>
6 : #include <algorithm>
7 : #include <vector>
8 :
9 : #include "spdlog/fmt/fmt.h"
10 :
11 : namespace elsa
12 : {
13 : using IndexRange = std::pair<int, int>;
14 :
15 : /*
16 : * @brief This is a iterable view over a Cartesian index space. This defines a rectangular
17 : * region of integer indices. This is mostly equivalent to a n-dimensional loop, e.g.:
18 : *
19 : * ```cpp
20 : * for(int i = istart; i < istop; ++i) {
21 : * for(int j = jstart; j < jstop; ++j) {
22 : * // nested arbitrarily deep
23 : * }
24 : * }
25 : * ```
26 : *
27 : * If you want to iterate over all indices form `0` to `upper` in a n-dimensional setting you
28 : * can:
29 : *
30 : * ```cpp
31 : * for(auto pos : CartesianIndex(upper)) {
32 : * // nested arbitrarily deep
33 : * }
34 : * ```
35 : *
36 : * Note, that this is as typically an exclusive range for the end, i.e. [lower, upper).
37 : * At no point any coefficient of `pos` take the value of `upper` (lower[i] <= pos[i] < upper[i]
38 : * for all i in 0, ... pos.size()).
39 : *
40 : * Assuming `x` is a vector of integers. If you don't want to start at `0`. Pass the lower bound
41 : * index as the first argument:
42 : *
43 : * ```cpp
44 : * for(auto pos : CartesianIndex(lower, upper)) {
45 : * // nested arbitrarily deep
46 : * }
47 : * ```
48 : *
49 : * The behaviour is quite similar to `std::views::iota` in C++20, but extended to n dimensions.
50 : * In this sense it is quite similar to Julies, check this
51 : * [blob](https://julialang.org/blog/2016/02/iteration/) for more information. This was also
52 : * the initial inspiration for this implementation.
53 : *
54 : */
55 : class CartesianIndices
56 : {
57 : private:
58 : index_t size_;
59 : IndexVector_t first_;
60 : IndexVector_t last_;
61 :
62 : // Tag struct, just to indicate the end iterator as sentinel. TODO: With a switch to C++
63 : // 20, remove this and let `end()` return a proper sentinel type
64 : struct as_sentinel {
65 : };
66 :
67 : static constexpr as_sentinel as_sentinel_t{};
68 :
69 : public:
70 : /// Vector of pairs, which represent the range for each dimension
71 : CartesianIndices(std::vector<IndexRange> ranges);
72 :
73 : /// Create an index space ranging from `0`, to `to` with the dimension of `to`
74 : CartesianIndices(const IndexVector_t& to);
75 :
76 : template <typename Vector>
77 : CartesianIndices(const Vector& to)
78 : : size_(to.size()), first_(IndexVector_t::Zero(size_)), last_(to.size())
79 31 : {
80 31 : using std::begin;
81 31 : std::copy_n(begin(to), size_, last_.begin());
82 31 : }
83 :
84 : /// Create an index space ranging from `from`, to `to` with the dimension of `ranges`
85 : CartesianIndices(const IndexVector_t& from, const IndexVector_t& to);
86 :
87 : template <typename Vector>
88 : CartesianIndices(const Vector& from, const Vector& to)
89 : : size_(to.size()), first_(from.size()), last_(to.size())
90 13 : {
91 13 : using std::begin;
92 13 : std::copy_n(begin(from), size_, first_.begin());
93 13 : std::copy_n(begin(to), size_, last_.begin());
94 13 : }
95 :
96 : /// Return the dimension of index space
97 : auto dims() const -> index_t;
98 :
99 : /// Return the number of coordinates of the space
100 : auto size() const -> index_t;
101 :
102 : /// return a range for a given dimension (0 <= i < dims())
103 : auto range(index_t i) const -> IndexRange;
104 :
105 : /// Return the lower bound of the index space
106 : auto first() -> IndexVector_t&;
107 : /// @overload
108 : auto first() const -> const IndexVector_t&;
109 :
110 : /// Return the upper bound of the index space
111 : auto last() -> IndexVector_t&;
112 : /// @overload
113 : auto last() const -> const IndexVector_t&;
114 :
115 : /// Random Access Iterator for index space
116 : struct iterator {
117 : private:
118 : IndexVector_t cur_;
119 : const IndexVector_t& begins_;
120 : const IndexVector_t& ends_;
121 :
122 : public:
123 : using value_type = IndexVector_t;
124 : using reference = value_type&;
125 : using pointer = value_type*;
126 : using difference_type = std::ptrdiff_t;
127 : using iterator_category = std::random_access_iterator_tag;
128 :
129 : iterator(as_sentinel, const IndexVector_t& begin, const IndexVector_t& end);
130 :
131 : explicit iterator(const IndexVector_t& vec);
132 :
133 : iterator(const IndexVector_t& cur, const IndexVector_t& begin,
134 : const IndexVector_t& end);
135 :
136 : // Dereference
137 : const IndexVector_t& operator*() const;
138 :
139 : // Increment
140 : iterator& operator++();
141 : iterator operator++(int);
142 :
143 : // decrement
144 : iterator& operator--();
145 : iterator operator--(int);
146 :
147 : // advance
148 : auto operator+=(difference_type n) -> CartesianIndices::iterator&;
149 : auto operator-=(difference_type n) -> CartesianIndices::iterator&;
150 :
151 : friend auto operator+(const iterator& iter, difference_type n) -> iterator;
152 : friend auto operator+(difference_type n, const iterator& iter) -> iterator;
153 :
154 : friend auto operator-(const iterator& iter, difference_type n) -> iterator;
155 : friend auto operator-(const iterator& lhs, const iterator& rhs) -> difference_type;
156 :
157 : auto operator[](difference_type n) const -> value_type;
158 :
159 : // comparison
160 : friend auto operator==(const iterator& lhs, const iterator& rhs) -> bool
161 208879 : {
162 208879 : return lhs.cur_ == rhs.cur_;
163 208879 : }
164 :
165 : friend auto operator!=(const iterator& lhs, const iterator& rhs) -> bool
166 208864 : {
167 208864 : return !(lhs == rhs);
168 208864 : }
169 :
170 : // relational operators
171 : friend auto operator<(const iterator& lhs, const iterator& rhs) -> bool;
172 : friend auto operator>(const iterator& lhs, const iterator& rhs) -> bool;
173 : friend auto operator<=(const iterator& lhs, const iterator& rhs) -> bool;
174 : friend auto operator>=(const iterator& lhs, const iterator& rhs) -> bool;
175 :
176 : private:
177 : auto at_end() -> bool;
178 : auto at_end() const -> bool;
179 :
180 : // TODO: Prefer a for loop, and remove recursive function
181 : auto inc_recursive(index_t N) -> void;
182 : auto inc() -> void;
183 :
184 : auto prev_recursive(index_t N) -> void;
185 : auto prev() -> void;
186 :
187 : auto distance_to_recusrive(const iterator& iter, difference_type N) const
188 : -> difference_type;
189 : auto distance_to(const iterator& iter) const -> difference_type;
190 :
191 : auto advance_recursive(difference_type n, index_t N) -> void;
192 : auto advance(difference_type n) -> void;
193 : };
194 :
195 : /// Return the begin iterator to the index space
196 : auto begin() -> iterator;
197 : /// Return the end iterator to the index space
198 : auto end() -> iterator;
199 :
200 : /// @overload
201 : auto begin() const -> iterator;
202 : /// @overload
203 : auto end() const -> iterator;
204 : };
205 :
206 : /// @brief Visit all neighbours with a certain distance `dist`, which have the same x-axis (i.e.
207 : /// they lie in the same y-z plane)
208 : ///
209 : /// The returned position can lie outside of the volume or have negative indices, so be careful
210 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist);
211 :
212 : /// @overload
213 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist);
214 :
215 : /// @brief Visit all neighbours with a certain distance `dist`, which have the same x-axis (i.e.
216 : /// they lie in the same y-z plane), plus make sure we stay in bounds of `lower`, and `upper`.
217 : ///
218 : /// Given a certain voxel in a volume, you can visit all close by voxels, which are still in the
219 : /// volume
220 : ///
221 : /// ```cpp
222 : /// auto volmin = {0, 0, 0};
223 : /// auto volmax = {128, 128, 128};
224 : /// auto current_pos = ...;
225 : /// for(auto neighbours : neighbours_in_slice(current_pos, 1, volmin, volmax)) {
226 : /// // ...
227 : /// }
228 : /// ```
229 : /// Note, that the current_pos is also visited!
230 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist,
231 : const IndexVector_t& lower, const IndexVector_t& upper);
232 :
233 : /// @overload
234 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist, index_t leadingDim,
235 : const IndexVector_t& lower, const IndexVector_t& upper);
236 : } // namespace elsa
237 :
238 : template <>
239 : struct fmt::formatter<elsa::CartesianIndices> {
240 : template <typename ParseContext>
241 : constexpr auto parse(ParseContext& ctx)
242 5 : {
243 5 : return ctx.begin();
244 5 : }
245 :
246 : template <typename FormatContext>
247 : auto format(const elsa::CartesianIndices& idx, FormatContext& ctx)
248 5 : {
249 5 : fmt::format_to(ctx.out(), "(");
250 :
251 9 : for (int i = 0; i < idx.dims() - 1; ++i) {
252 4 : auto p = idx.range(i);
253 4 : fmt::format_to(ctx.out(), "{}:{}, ", p.first, p.second);
254 4 : }
255 5 : auto p = idx.range(idx.dims() - 1);
256 5 : return fmt::format_to(ctx.out(), "{}:{})", p.first, p.second);
257 5 : }
258 : };
|
