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