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