1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
/* Copyright (c) 2008-2022 the MRtrix3 contributors.
*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
* Covered Software is provided under this License on an "as is"
* basis, without warranty of any kind, either expressed, implied, or
* statutory, including, without limitation, warranties that the
* Covered Software is free of defects, merchantable, fit for a
* particular purpose or non-infringing.
* See the Mozilla Public License v. 2.0 for more details.
*
* For more details, see http://www.mrtrix.org/.
*/
#ifndef __algo_neighbourhooditerator_h__
#define __algo_neighbourhooditerator_h__
#include "types.h"
#include "algo/iterator.h"
namespace MR
{
/** \defgroup loop Looping functions
@{ */
//! a dummy image to iterate over a certain neighbourhood, useful for multi-threaded looping.
// Does not work properly with Loop() functions! Use instead:
//
// vector extent(iter.ndim(),3) // max number of voxels to iterate over
// auto niter = NeighbourhoodIterator(iter, extent);
// while(niter.loop()){
// std::cerr << niter << std::endl;
// }
//
//
class NeighbourhoodIterator
{ MEMALIGN(NeighbourhoodIterator)
public:
NeighbourhoodIterator() = delete;
template <class IteratorType>
NeighbourhoodIterator (const IteratorType& iter, const vector<size_t>& extent) :
dim (iter.ndim()),
offset (iter.ndim()),
// pos (iter.ndim()),
pos_orig (iter.ndim()),
ext (container_cast<decltype(ext)>(extent)),
has_next_(false) {
assert (iter.ndim() == extent.size());
pos.resize(iter.ndim());
for (size_t i = 0; i < iter.ndim(); ++i){
ext[i] = (ext[i]-1) / 2;
offset[i] = iter.index(i);
// set pos to lower bound
pos[i] = (offset[i] - ext[i] < 0) ? 0 : offset[i] - ext[i];
pos_orig[i] = pos[i];
// upper bound:
auto high = (offset[i] + ext[i] >= iter.size(i)) ? iter.size(i) - 1 : offset[i] + ext[i];
dim[i] = high - pos[i] + 1;
}
}
size_t ndim () const { return dim.size(); }
ssize_t size (size_t axis) const { return dim[axis]; }
const ssize_t& index (size_t axis) const { return pos[axis]; }
ssize_t& index (size_t axis) { return pos[axis]; }
const Eigen::Matrix< ssize_t, 1, Eigen::Dynamic > get_pos() const { return pos; }
const ssize_t& extent (size_t axis) const { return dim[axis]; }
const ssize_t& centre (size_t axis) const { return offset[axis]; }
void reset (size_t axis) { pos[axis] = pos_orig[axis]; }
bool loop () {
if (not this->has_next_){
this->has_next_ = true;
for (auto axis = dim.size(); axis-- != 0; )
this->reset(axis);
return true;
}
for (auto axis = dim.size(); axis-- != 0; ) { // TODO loop in stride order
++pos[axis];
if (pos[axis] != (pos_orig[axis] + dim[axis]))
return true;
this->reset(axis);
}
this->has_next_ = false;
return false;
}
friend std::ostream& operator<< (std::ostream& stream, const NeighbourhoodIterator& V) {
stream << "neighbourhood iterator, position [ ";
for (size_t n = 0; n < V.ndim(); ++n)
stream << V.index(n) << " ";
stream << "], extent [ ";
for (size_t n = 0; n < V.ndim(); ++n)
stream << V.extent(n) << " ";
stream << "], centre [ ";
for (size_t n = 0; n < V.ndim(); ++n)
stream << V.centre(n) << " ";
stream << "]";
return stream;
}
private:
vector<ssize_t> dim, offset, pos_orig, ext;
Eigen::Matrix< ssize_t, 1, Eigen::Dynamic > pos;
bool has_next_;
void value () const { assert (0); }
};
//! @}
}
#endif
|
