1 #ifndef MAHOTAS_FILTER_H_INCLUDE_GUARD_
2 #define MAHOTAS_FILTER_H_INCLUDE_GUARD_
3 // Copyright (C) 2003-2005 Peter J. Verveer
4 // Copyright (C) 2010-2013 Luis Pedro Coelho
5 // LICENSE: MIT
6
7 #include <vector>
8 #include <cassert>
9 #include <limits>
10 #include "numpypp/array.hpp"
11
12
13 /* The different boundary conditions. The mirror condition is not used
14 by the python code, but C code is kept around in case we might wish
15 to add it. */
16 typedef enum {
17 ExtendNearest = 0,
18 ExtendWrap = 1,
19 ExtendReflect = 2,
20 ExtendMirror = 3,
21 ExtendConstant = 4,
22 ExtendIgnore = 5
23 } ExtendMode;
24 const ExtendMode ExtendLast = ExtendIgnore;
25
26 const npy_intp border_flag_value = std::numeric_limits<npy_intp>::max();
27
28 inline
fix_offset(const ExtendMode mode,npy_intp cc,const npy_intp len)29 npy_intp fix_offset(const ExtendMode mode, npy_intp cc, const npy_intp len) {
30 /* apply boundary conditions, if necessary: */
31 switch (mode) {
32 case ExtendMirror:
33 if (cc < 0) {
34 if (len <= 1) {
35 return 0;
36 } else {
37 int sz2 = 2 * len - 2;
38 cc = sz2 * (int)(-cc / sz2) + cc;
39 return cc <= 1 - len ? cc + sz2 : -cc;
40 }
41 } else if (cc >= len) {
42 if (len <= 1) {
43 return 0;
44 } else {
45 int sz2 = 2 * len - 2;
46 cc -= sz2 * (int)(cc / sz2);
47 if (cc >= len)
48 cc = sz2 - cc;
49 }
50 }
51 return cc;
52
53 case ExtendReflect:
54 if (cc < 0) {
55 if (len <= 1) {
56 return 0;
57 } else {
58 int sz2 = 2 * len;
59 if (cc < -sz2)
60 cc = sz2 * (int)(-cc / sz2) + cc;
61 cc = cc < -len ? cc + sz2 : -cc - 1;
62 }
63 } else if (cc >= len) {
64 if (len <= 1) {
65 return 0;
66 } else {
67 int sz2 = 2 * len;
68 cc -= sz2 * (int)(cc / sz2);
69 if (cc >= len)
70 cc = sz2 - cc - 1;
71 }
72 }
73 return cc;
74 case ExtendWrap:
75 if (cc < 0) {
76 if (len <= 1) {
77 return 0;
78 } else {
79 int sz = len;
80 cc += sz * (int)(-cc / sz);
81 if (cc < 0)
82 cc += sz;
83 }
84 } else if (cc >= len) {
85 if (len <= 1) {
86 return 0;
87 } else {
88 int sz = len;
89 cc -= sz * (int)(cc / sz);
90 }
91 }
92 return cc;
93 case ExtendNearest:
94 if (cc < 0) return 0;
95 if (cc >= len) return len - 1;
96 return cc;
97 case ExtendIgnore:
98 case ExtendConstant:
99 if (cc < 0 || cc >= len)
100 return border_flag_value;
101 return cc;
102 }
103 assert(false); // We should never get here
104 return 0;
105 }
106
107
108 int init_filter_offsets(PyArrayObject *array, bool *footprint,
109 const npy_intp * const fshape, npy_intp* origins,
110 const ExtendMode mode, std::vector<npy_intp>& offsets,
111 std::vector<npy_intp>* coordinate_offsets);
112 void init_filter_iterator(const int rank, const npy_intp *fshape,
113 const npy_intp filter_size, const npy_intp *ashape,
114 const npy_intp *origins,
115 npy_intp* strides, npy_intp* backstrides,
116 npy_intp* minbound, npy_intp* maxbound);
117
118 template <typename T>
119 struct filter_iterator {
120 /* Move to the next point in an array, possible changing the pointer
121 to the filter offsets when moving into a different region in the
122 array: */
123 filter_iterator(PyArrayObject* array, PyArrayObject* filter, ExtendMode mode=ExtendNearest, bool compress=true)
filter_data_filter_iterator124 :filter_data_(numpy::ndarray_cast<T*>(filter))
125 ,own_filter_data_(false)
126 ,nd_(PyArray_NDIM(array))
127 {
128 numpy::aligned_array<T> filter_array(filter);
129 const npy_intp filter_size = filter_array.size();
130 bool* footprint = 0;
131 if (compress) {
132 footprint = new bool[filter_size];
133 typename numpy::aligned_array<T>::iterator fiter = filter_array.begin();
134 for (int i = 0; i != filter_size; ++i, ++fiter) {
135 footprint[i] = !!(*fiter);
136 }
137 }
138 size_ = init_filter_offsets(array, footprint, PyArray_DIMS(filter), 0,
139 mode, offsets_, 0);
140 if (compress) {
141 int j = 0;
142 T* new_filter_data = new T[size_];
143 typename numpy::aligned_array<T>::iterator fiter = filter_array.begin();
144 for (int i = 0; i != filter_size; ++i, ++fiter) {
145 if (*fiter) {
146 new_filter_data[j++] = *fiter;
147 }
148 }
149 filter_data_ = new_filter_data;
150 own_filter_data_ = true;
151 delete [] footprint;
152 }
153
154 init_filter_iterator(PyArray_NDIM(filter), PyArray_DIMS(filter), size_,
155 PyArray_DIMS(array), /*origins*/0,
156 this->strides_, this->backstrides_,
157 this->minbound_, this->maxbound_);
158 cur_offsets_idx_ = this->offsets_.begin();
159 }
~filter_iteratorfilter_iterator160 ~filter_iterator() {
161 if (own_filter_data_) delete [] filter_data_;
162 }
163 template <typename OtherIterator>
iterate_bothfilter_iterator164 void iterate_both(OtherIterator& iterator) {
165 for (int d = 0; d < nd_; ++d) {
166 const npy_intp p = iterator.index_rev(d);
167 if (p < (iterator.dimension_rev(d) - 1)) {
168 if (p < this->minbound_[d] || p >= this->maxbound_[d]) {
169 this->cur_offsets_idx_ += this->strides_[d];
170 }
171 break;
172 }
173 this->cur_offsets_idx_ -= this->backstrides_[d];
174 assert(this->cur_offsets_idx_ >= this->offsets_.begin());
175 assert(this->cur_offsets_idx_ < this->offsets_.end());
176 }
177 ++iterator;
178 }
179
180 template <typename OtherIterator>
retrievefilter_iterator181 bool retrieve(const OtherIterator& iterator, const npy_intp j, T& array_val) {
182 assert((j >= 0) && (j < size_));
183 if (this->cur_offsets_idx_[j] == border_flag_value) return false;
184 array_val = *( (&*iterator) + this->cur_offsets_idx_[j]);
185 return true;
186 }
187 template <typename OtherIterator>
setfilter_iterator188 void set(const OtherIterator& iterator, npy_intp j, const T& val) {
189 assert(this->cur_offsets_idx_[j] != border_flag_value);
190 *( (&*iterator) + this->cur_offsets_idx_[j]) = val;
191 }
192
193 const T& operator [] (const npy_intp j) const { assert(j < size_); return filter_data_[j]; }
sizefilter_iterator194 npy_intp size() const { return size_; }
195 private:
196 const T* filter_data_;
197 bool own_filter_data_;
198 std::vector<npy_intp>::const_iterator cur_offsets_idx_;
199 npy_intp size_;
200 const npy_intp nd_;
201 std::vector<npy_intp> offsets_;
202 npy_intp strides_[NPY_MAXDIMS];
203 npy_intp backstrides_[NPY_MAXDIMS];
204 npy_intp minbound_[NPY_MAXDIMS];
205 npy_intp maxbound_[NPY_MAXDIMS];
206 };
207
208 #endif // MAHOTAS_FILTER_H_INCLUDE_GUARD_
209