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 0 : CartesianIndices::CartesianIndices(std::vector<IndexRange> ranges) : size_(ranges.size()) {}
26 :
27 0 : CartesianIndices::CartesianIndices(const IndexVector_t& to)
28 0 : : size_(to.size()), first_(IndexVector_t::Zero(size_)), last_(to)
29 : {
30 0 : }
31 :
32 0 : CartesianIndices::CartesianIndices(const IndexVector_t& from, const IndexVector_t& to)
33 0 : : size_(to.size()), first_(from), last_(to)
34 : {
35 0 : if (from.size() != to.size()) {
36 0 : throw InvalidArgumentError("CartesianIndices: vectors must be of same size");
37 : }
38 0 : }
39 :
40 0 : auto CartesianIndices::dims() const -> index_t { return size_; }
41 :
42 0 : auto CartesianIndices::size() const -> index_t { return (last_ - first_).prod(); }
43 :
44 0 : auto CartesianIndices::range(index_t i) const -> IndexRange { return {first_[i], last_[i]}; }
45 :
46 0 : auto CartesianIndices::first() -> IndexVector_t& { return first_; }
47 :
48 0 : auto CartesianIndices::first() const -> const IndexVector_t& { return first_; }
49 :
50 0 : auto CartesianIndices::last() -> IndexVector_t& { return last_; }
51 :
52 0 : auto CartesianIndices::last() const -> const IndexVector_t& { return last_; }
53 :
54 0 : auto CartesianIndices::begin() -> iterator { return iterator(first(), first(), last()); }
55 :
56 0 : auto CartesianIndices::begin() const -> iterator { return {first(), first(), last()}; }
57 :
58 0 : auto CartesianIndices::end() -> iterator { return {as_sentinel_t, first(), last()}; }
59 :
60 0 : auto CartesianIndices::end() const -> iterator { return {as_sentinel_t, first(), last()}; }
61 :
62 : /****** Iterator implementation ******/
63 0 : CartesianIndices::iterator::iterator(as_sentinel, const IndexVector_t& begin,
64 0 : const IndexVector_t& end)
65 0 : : cur_(end), begins_(end), ends_(end)
66 : {
67 0 : cur_.tail(cur_.size() - 1) = begin.tail(cur_.size() - 1);
68 0 : }
69 :
70 0 : CartesianIndices::iterator::iterator(const IndexVector_t& vec)
71 0 : : cur_(vec), begins_(vec), ends_(vec)
72 : {
73 0 : }
74 :
75 0 : CartesianIndices::iterator::iterator(const IndexVector_t& cur, const IndexVector_t& begin,
76 0 : const IndexVector_t& end)
77 0 : : cur_(cur), begins_(begin), ends_(end)
78 : {
79 0 : }
80 :
81 0 : auto CartesianIndices::iterator::operator*() const -> const IndexVector_t& { return cur_; }
82 :
83 0 : auto CartesianIndices::iterator::operator++() -> iterator&
84 : {
85 0 : inc();
86 0 : return *this;
87 : }
88 :
89 0 : auto CartesianIndices::iterator::operator++(int) -> iterator
90 : {
91 0 : auto copy = *this;
92 0 : ++(*this);
93 0 : return copy;
94 0 : }
95 :
96 0 : auto CartesianIndices::iterator::operator--() -> iterator&
97 : {
98 0 : prev();
99 0 : return *this;
100 : }
101 :
102 0 : auto CartesianIndices::iterator::operator--(int) -> iterator
103 : {
104 0 : auto copy = *this;
105 0 : --(*this);
106 0 : return copy;
107 0 : }
108 :
109 0 : auto CartesianIndices::iterator::operator+=(difference_type n) -> iterator&
110 : {
111 0 : advance(n);
112 0 : return *this;
113 : }
114 :
115 0 : auto CartesianIndices::iterator::operator-=(difference_type n) -> iterator&
116 : {
117 0 : advance(-n);
118 0 : return *this;
119 : }
120 :
121 0 : auto operator+(const CartesianIndices::iterator& iter,
122 : CartesianIndices::iterator::difference_type n) -> CartesianIndices::iterator
123 : {
124 0 : auto copy = iter;
125 0 : copy += n;
126 0 : return copy;
127 0 : }
128 :
129 0 : auto operator+(CartesianIndices::iterator::difference_type n,
130 : const CartesianIndices::iterator& iter) -> CartesianIndices::iterator
131 : {
132 0 : return iter + n;
133 : }
134 :
135 0 : auto operator-(const CartesianIndices::iterator& iter,
136 : CartesianIndices::iterator::difference_type n) -> CartesianIndices::iterator
137 : {
138 0 : auto copy = iter;
139 0 : copy -= n;
140 0 : return copy;
141 0 : }
142 :
143 0 : auto operator-(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
144 : -> CartesianIndices::iterator::difference_type
145 : {
146 0 : return rhs.distance_to(lhs);
147 : }
148 :
149 0 : auto CartesianIndices::iterator::operator[](difference_type n) const -> value_type
150 : {
151 : // TODO: Make this more efficient
152 0 : auto copy = *this;
153 0 : copy += n;
154 0 : return *copy;
155 0 : }
156 :
157 0 : auto operator<(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
158 : -> bool
159 : {
160 : // Not sure why I have to instantiate them this way, else it doesn't work...
161 : using Array = Eigen::Array<index_t, Eigen::Dynamic, 1>;
162 0 : const Array a1 = lhs.cur_.array();
163 0 : const Array a2 = rhs.cur_.array();
164 :
165 0 : return (a1 < a2).any();
166 0 : }
167 :
168 0 : auto operator>(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
169 : -> bool
170 : {
171 0 : return rhs < lhs;
172 : }
173 :
174 0 : auto operator<=(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
175 : -> bool
176 : {
177 0 : return !(lhs > rhs);
178 : }
179 :
180 0 : auto operator>=(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
181 : -> bool
182 : {
183 0 : return !(lhs < rhs);
184 : }
185 :
186 0 : auto CartesianIndices::iterator::at_end() -> bool { return cur_[0] == ends_[0]; }
187 :
188 0 : auto CartesianIndices::iterator::at_end() const -> bool { return cur_[0] == ends_[0]; }
189 :
190 0 : auto CartesianIndices::iterator::inc_recursive(index_t N) -> void
191 : {
192 0 : auto& iter = cur_[N];
193 :
194 0 : if (++iter == ends_[N] && N > 0) {
195 0 : iter = begins_[N];
196 0 : inc_recursive(N - 1);
197 : }
198 0 : }
199 :
200 0 : auto CartesianIndices::iterator::inc() -> void { inc_recursive(cur_.size() - 1); }
201 :
202 0 : auto CartesianIndices::iterator::prev_recursive(index_t N) -> void
203 : {
204 0 : auto& iter = cur_[N];
205 :
206 0 : if (iter == begins_[N] && N > 0) {
207 0 : iter = ends_[N];
208 0 : inc_recursive(N - 1);
209 : }
210 0 : --iter;
211 0 : }
212 :
213 0 : auto CartesianIndices::iterator::prev() -> void { prev_recursive(cur_.size() - 1); }
214 :
215 0 : auto CartesianIndices::iterator::distance_to_recusrive(const iterator& other,
216 : difference_type N) const
217 : -> difference_type
218 : {
219 0 : if (N == 0) {
220 0 : return other.cur_[N] - cur_[N];
221 : } else {
222 0 : const auto d = distance_to_recusrive(other, N - 1);
223 0 : const auto scale = ends_[N] - begins_[N];
224 0 : const auto increment = other.cur_[N] - cur_[N];
225 :
226 0 : return difference_type{(d * scale) + increment};
227 : }
228 : }
229 :
230 0 : auto CartesianIndices::iterator::distance_to(const iterator& other) const -> difference_type
231 : {
232 0 : auto N = cur_.size() - 1;
233 0 : return distance_to_recusrive(other, N);
234 : }
235 :
236 0 : auto CartesianIndices::iterator::advance_recursive(difference_type n, index_t N) -> void
237 : {
238 0 : if (n == 0) {
239 0 : return;
240 : }
241 :
242 0 : auto& iter = cur_[N];
243 0 : const auto size = ends_[N] - begins_[N];
244 0 : const auto first = begins_[N];
245 :
246 0 : auto const idx = as<difference_type>(iter - first);
247 0 : n += idx;
248 :
249 0 : auto div = size ? n / size : 0;
250 0 : auto mod = size ? n % size : 0;
251 :
252 0 : if (N != 0) {
253 0 : if (mod < 0) {
254 0 : mod += size;
255 0 : div--;
256 : }
257 0 : advance_recursive(div, N - 1);
258 : } else {
259 0 : if (div > 0) {
260 0 : mod = size;
261 : }
262 : }
263 :
264 0 : iter = first + mod;
265 : }
266 :
267 0 : auto CartesianIndices::iterator::advance(difference_type n) -> void
268 : {
269 0 : advance_recursive(n, cur_.size() - 1);
270 0 : }
271 :
272 0 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist)
273 : {
274 0 : Eigen::IOFormat format(4, 0, ", ", "", "", "", "[", "]");
275 0 : IndexVector_t from = pos;
276 0 : from.tail(pos.size() - 1).array() -= dist.array();
277 :
278 0 : IndexVector_t to = pos;
279 0 : to.tail(pos.size() - 1).array() += dist.array() + 1;
280 0 : to[0] += 1; // Be sure this is incremented, so we actually iterate over it
281 :
282 0 : return {from, to};
283 0 : }
284 :
285 0 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist)
286 : {
287 0 : return neighbours_in_slice(pos, IndexVector_t::Constant(pos.size() - 1, dist));
288 : }
289 :
290 0 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist,
291 : const IndexVector_t& lower, const IndexVector_t& upper)
292 : {
293 0 : IndexVector_t from = pos;
294 0 : from.array() -= dist.array();
295 0 : from = from.cwiseMax(lower);
296 :
297 0 : IndexVector_t to = pos;
298 0 : to.array() += dist.array() + 1;
299 0 : to = to.cwiseMin(upper);
300 :
301 : // FIXME: sometimes this happens when padding is added to an aabb
302 0 : from = (to.array() == from.array()).select(from.array() - 1, from);
303 :
304 0 : return {from, to};
305 0 : }
306 :
307 0 : CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist, index_t leadingDim,
308 : const IndexVector_t& lower, const IndexVector_t& upper)
309 : {
310 0 : IndexVector_t distvec = IndexVector_t::Constant(pos.size(), dist);
311 0 : distvec[leadingDim] = 0;
312 :
313 0 : return neighbours_in_slice(pos, distvec, lower, upper);
314 0 : }
315 : } // namespace elsa
|
