Line data Source code
1 : #include "CartesianIndices.h"
2 : #include "TypeCasts.hpp"
3 : #include <iostream>
4 : #include <iterator>
5 : #include <numeric>
6 :
7 : namespace elsa
8 : {
9 : template <typename Fn>
10 : auto map(std::vector<IndexRange> v, Fn fn)
11 : {
12 : IndexVector_t result(v.size());
13 : std::transform(v.begin(), v.end(), result.begin(), std::move(fn));
14 : return result;
15 : }
16 :
17 : template <std::size_t N, typename Fn>
18 : auto map(std::array<IndexRange, N> v, index_t size, Fn fn)
19 : {
20 : IndexVector_t result(size);
21 : std::transform(v.begin(), v.begin() + size, result.begin(), std::move(fn));
22 : return result;
23 : }
24 :
25 : CartesianIndices::CartesianIndices(std::vector<IndexRange> ranges)
26 : : size_(static_cast<index_t>(ranges.size()))
27 0 : {
28 0 : }
29 :
30 : CartesianIndices::CartesianIndices(const IndexVector_t& to)
31 : : size_(to.size()), first_(IndexVector_t::Zero(size_)), last_(to)
32 1663 : {
33 1663 : }
34 :
35 : CartesianIndices::CartesianIndices(const IndexVector_t& from, const IndexVector_t& to)
36 : : size_(to.size()), first_(from), last_(to)
37 7970 : {
38 7970 : if (from.size() != to.size()) {
39 0 : throw InvalidArgumentError("CartesianIndices: vectors must be of same size");
40 0 : }
41 7970 : }
42 :
43 : auto CartesianIndices::dims() const -> index_t
44 20 : {
45 20 : return size_;
46 20 : }
47 :
48 : auto CartesianIndices::size() const -> index_t
49 491 : {
50 491 : return (last_ - first_).prod();
51 491 : }
52 :
53 : auto CartesianIndices::range(index_t i) const -> IndexRange
54 9 : {
55 9 : return {first_[i], last_[i]};
56 9 : }
57 :
58 : auto CartesianIndices::first() -> IndexVector_t&
59 24088 : {
60 24088 : return first_;
61 24088 : }
62 :
63 : auto CartesianIndices::first() const -> const IndexVector_t&
64 66 : {
65 66 : return first_;
66 66 : }
67 :
68 : auto CartesianIndices::last() -> IndexVector_t&
69 16061 : {
70 16061 : return last_;
71 16061 : }
72 :
73 : auto CartesianIndices::last() const -> const IndexVector_t&
74 40 : {
75 40 : return last_;
76 40 : }
77 :
78 : auto CartesianIndices::begin() -> iterator
79 8039 : {
80 8039 : return iterator(first(), first(), last());
81 8039 : }
82 :
83 : auto CartesianIndices::begin() const -> iterator
84 26 : {
85 26 : return {first(), first(), last()};
86 26 : }
87 :
88 : auto CartesianIndices::end() -> iterator
89 8033 : {
90 8033 : return {as_sentinel_t, first(), last()};
91 8033 : }
92 :
93 : auto CartesianIndices::end() const -> iterator
94 11 : {
95 11 : return {as_sentinel_t, first(), last()};
96 11 : }
97 :
98 : /****** Iterator implementation ******/
99 : CartesianIndices::iterator::iterator(as_sentinel, const IndexVector_t& begin,
100 : const IndexVector_t& end)
101 : : cur_(end), begins_(end), ends_(end)
102 8026 : {
103 8026 : cur_.tail(cur_.size() - 1) = begin.tail(cur_.size() - 1);
104 8026 : }
105 :
106 : CartesianIndices::iterator::iterator(const IndexVector_t& vec)
107 : : cur_(vec), begins_(vec), ends_(vec)
108 0 : {
109 0 : }
110 :
111 : CartesianIndices::iterator::iterator(const IndexVector_t& cur, const IndexVector_t& begin,
112 : const IndexVector_t& end)
113 : : cur_(cur), begins_(begin), ends_(end)
114 8061 : {
115 8061 : }
116 :
117 : auto CartesianIndices::iterator::operator*() const -> const IndexVector_t&
118 200941 : {
119 200941 : return cur_;
120 200941 : }
121 :
122 : auto CartesianIndices::iterator::operator++() -> iterator&
123 201389 : {
124 201389 : inc();
125 201389 : return *this;
126 201389 : }
127 :
128 : auto CartesianIndices::iterator::operator++(int) -> iterator
129 1 : {
130 1 : auto copy = *this;
131 1 : ++(*this);
132 1 : return copy;
133 1 : }
134 :
135 : auto CartesianIndices::iterator::operator--() -> iterator&
136 2 : {
137 2 : prev();
138 2 : return *this;
139 2 : }
140 :
141 : auto CartesianIndices::iterator::operator--(int) -> iterator
142 1 : {
143 1 : auto copy = *this;
144 1 : --(*this);
145 1 : return copy;
146 1 : }
147 :
148 : auto CartesianIndices::iterator::operator+=(difference_type n) -> iterator&
149 34 : {
150 34 : advance(n);
151 34 : return *this;
152 34 : }
153 :
154 : auto CartesianIndices::iterator::operator-=(difference_type n) -> iterator&
155 14 : {
156 14 : advance(-n);
157 14 : return *this;
158 14 : }
159 :
160 : auto operator+(const CartesianIndices::iterator& iter,
161 : CartesianIndices::iterator::difference_type n) -> CartesianIndices::iterator
162 12 : {
163 12 : auto copy = iter;
164 12 : copy += n;
165 12 : return copy;
166 12 : }
167 :
168 : auto operator+(CartesianIndices::iterator::difference_type n,
169 : const CartesianIndices::iterator& iter) -> CartesianIndices::iterator
170 6 : {
171 6 : return iter + n;
172 6 : }
173 :
174 : auto operator-(const CartesianIndices::iterator& iter,
175 : CartesianIndices::iterator::difference_type n) -> CartesianIndices::iterator
176 6 : {
177 6 : auto copy = iter;
178 6 : copy -= n;
179 6 : return copy;
180 6 : }
181 :
182 : auto operator-(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
183 : -> CartesianIndices::iterator::difference_type
184 35 : {
185 35 : return rhs.distance_to(lhs);
186 35 : }
187 :
188 : auto CartesianIndices::iterator::operator[](difference_type n) const -> value_type
189 5 : {
190 : // TODO: Make this more efficient
191 5 : auto copy = *this;
192 5 : copy += n;
193 5 : return *copy;
194 5 : }
195 :
196 : auto operator<(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
197 : -> bool
198 16 : {
199 : // Not sure why I have to instantiate them this way, else it doesn't work...
200 16 : using Array = Eigen::Array<index_t, Eigen::Dynamic, 1>;
201 16 : const Array a1 = lhs.cur_.array();
202 16 : const Array a2 = rhs.cur_.array();
203 :
204 16 : return (a1 < a2).any();
205 16 : }
206 :
207 : auto operator>(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
208 : -> bool
209 8 : {
210 8 : return rhs < lhs;
211 8 : }
212 :
213 : auto operator<=(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
214 : -> bool
215 4 : {
216 4 : return !(lhs > rhs);
217 4 : }
218 :
219 : auto operator>=(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
220 : -> bool
221 4 : {
222 4 : return !(lhs < rhs);
223 4 : }
224 :
225 : auto CartesianIndices::iterator::at_end() -> bool
226 0 : {
227 0 : return cur_[0] == ends_[0];
228 0 : }
229 :
230 : auto CartesianIndices::iterator::at_end() const -> bool
231 0 : {
232 0 : return cur_[0] == ends_[0];
233 0 : }
234 :
235 : auto CartesianIndices::iterator::inc_recursive(index_t N) -> void
236 244962 : {
237 244962 : auto& iter = cur_[N];
238 :
239 244962 : if (++iter == ends_[N] && N > 0) {
240 43610 : iter = begins_[N];
241 43610 : inc_recursive(N - 1);
242 43610 : }
243 244962 : }
244 :
245 : auto CartesianIndices::iterator::inc() -> void
246 201394 : {
247 201394 : inc_recursive(cur_.size() - 1);
248 201394 : }
249 :
250 : auto CartesianIndices::iterator::prev_recursive(index_t N) -> void
251 2 : {
252 2 : auto& iter = cur_[N];
253 :
254 2 : if (iter == begins_[N] && N > 0) {
255 0 : iter = ends_[N];
256 0 : inc_recursive(N - 1);
257 0 : }
258 2 : --iter;
259 2 : }
260 :
261 : auto CartesianIndices::iterator::prev() -> void
262 2 : {
263 2 : prev_recursive(cur_.size() - 1);
264 2 : }
265 :
266 : auto CartesianIndices::iterator::distance_to_recusrive(const iterator& other,
267 : difference_type N) const
268 : -> difference_type
269 70 : {
270 70 : if (N == 0) {
271 35 : return other.cur_[N] - cur_[N];
272 35 : } else {
273 35 : const auto d = distance_to_recusrive(other, N - 1);
274 35 : const auto scale = ends_[N] - begins_[N];
275 35 : const auto increment = other.cur_[N] - cur_[N];
276 :
277 35 : return difference_type{(d * scale) + increment};
278 35 : }
279 70 : }
280 :
281 : auto CartesianIndices::iterator::distance_to(const iterator& other) const -> difference_type
282 35 : {
283 35 : auto N = cur_.size() - 1;
284 35 : return distance_to_recusrive(other, N);
285 35 : }
286 :
287 : auto CartesianIndices::iterator::advance_recursive(difference_type n, index_t N) -> void
288 76 : {
289 76 : if (n == 0) {
290 18 : return;
291 18 : }
292 :
293 58 : auto& iter = cur_[N];
294 58 : const auto size = ends_[N] - begins_[N];
295 58 : const auto first = begins_[N];
296 :
297 58 : auto const idx = as<difference_type>(iter - first);
298 58 : n += idx;
299 :
300 58 : auto div = size ? n / size : 0;
301 58 : auto mod = size ? n % size : 0;
302 :
303 58 : if (N != 0) {
304 28 : if (mod < 0) {
305 0 : mod += size;
306 0 : div--;
307 0 : }
308 28 : advance_recursive(div, N - 1);
309 30 : } else {
310 30 : if (div > 0) {
311 7 : mod = size;
312 7 : }
313 30 : }
314 :
315 58 : iter = first + mod;
316 58 : }
317 :
318 : auto CartesianIndices::iterator::advance(difference_type n) -> void
319 48 : {
320 48 : advance_recursive(n, cur_.size() - 1);
321 48 : }
322 :
323 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist)
324 2 : {
325 2 : Eigen::IOFormat format(4, 0, ", ", "", "", "", "[", "]");
326 2 : IndexVector_t from = pos;
327 2 : from.tail(pos.size() - 1).array() -= dist.array();
328 :
329 2 : IndexVector_t to = pos;
330 2 : to.tail(pos.size() - 1).array() += dist.array() + 1;
331 2 : to[0] += 1; // Be sure this is incremented, so we actually iterate over it
332 :
333 2 : return {from, to};
334 2 : }
335 :
336 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist)
337 2 : {
338 2 : return neighbours_in_slice(pos, IndexVector_t::Constant(pos.size() - 1, dist));
339 2 : }
340 :
341 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist,
342 : const IndexVector_t& lower, const IndexVector_t& upper)
343 7951 : {
344 7951 : IndexVector_t from = pos;
345 7951 : from.array() -= dist.array();
346 7951 : from = from.cwiseMax(lower);
347 :
348 7951 : IndexVector_t to = pos;
349 7951 : to.array() += dist.array() + 1;
350 7951 : to = to.cwiseMin(upper);
351 :
352 : // FIXME: sometimes this happens when padding is added to an aabb
353 7951 : from = (to.array() == from.array()).select(from.array() - 1, from);
354 :
355 7951 : return {from, to};
356 7951 : }
357 :
358 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist, index_t leadingDim,
359 : const IndexVector_t& lower, const IndexVector_t& upper)
360 2 : {
361 2 : IndexVector_t distvec = IndexVector_t::Constant(pos.size(), dist);
362 2 : distvec[leadingDim] = 0;
363 :
364 2 : return neighbours_in_slice(pos, distvec, lower, upper);
365 2 : }
366 : } // namespace elsa
|
