1 // Copyright (C) 2008-2019, Luis Pedro Coelho <luis@luispedro.org>
2 // vim: set ts=4 sts=4 sw=4 expandtab smartindent:
3 //
4 // License: MIT
5 
6 #include <algorithm>
7 #include <queue>
8 #include <vector>
9 #include <cstdio>
10 #include <limits>
11 #include <iostream>
12 
13 #include "numpypp/array.hpp"
14 #include "numpypp/dispatch.hpp"
15 #include "utils.hpp"
16 
17 #include "_filters.h"
18 
19 namespace{
20 
21 const char TypeErrorMsg[] =
22     "Type not understood. "
23     "This is caused by either a direct call to _morph (which is dangerous: types are not checked!) or a bug in mahotas.\n";
24 
25 // Using std::abs fails on Windows as it is not defined for all integer types:
26 template <typename T>
t_abs(T val)27 inline T t_abs(T val) { return (val >= 0 ? val : -val); }
28 
29 template<typename T>
subm(numpy::aligned_array<T> a,const numpy::aligned_array<T> b)30 void subm(numpy::aligned_array<T> a, const numpy::aligned_array<T> b) {
31     gil_release nogil;
32     const numpy::index_type N = a.size();
33     typename numpy::aligned_array<T>::iterator ita = a.begin();
34     typename numpy::aligned_array<T>::const_iterator itb = b.begin();
35     for (numpy::index_type i = 0; i != N; ++i, ++ita, ++itb) {
36         if (std::numeric_limits<T>::is_signed) {
37             T val = *ita - *itb;
38             if (*itb >= 0 && (val <= *ita)) *ita = val; // subtracted a positive number, no underflow
39             else if (*itb < 0 && val > *ita) *ita = val; // subtracted a negative number, no overlow
40             else if (*itb >= 0) *ita = std::numeric_limits<T>::min(); // this is the underflow case
41             else *ita = std::numeric_limits<T>::max(); // this is the overflow case
42         } else {
43             if (*itb > *ita) *ita = T();
44             else *ita -= *itb;
45         }
46     }
47 }
48 
49 
py_subm(PyObject * self,PyObject * args)50 PyObject* py_subm(PyObject* self, PyObject* args) {
51     PyArrayObject* a;
52     PyArrayObject* b;
53     if (!PyArg_ParseTuple(args, "OO", &a, &b)) return NULL;
54     if (!numpy::are_arrays(a,b) ||
55         !numpy::same_shape(a, b) ||
56         !numpy::equiv_typenums(a,b)) {
57         PyErr_SetString(PyExc_RuntimeError, TypeErrorMsg);
58         return NULL;
59     }
60 #define HANDLE(type) \
61     subm<type>(numpy::aligned_array<type>(a), numpy::aligned_array<type>(b));
62     SAFE_SWITCH_ON_INTEGER_TYPES_OF(a);
63 #undef HANDLE
64 
65     Py_XINCREF(a);
66     return PyArray_Return(a);
67 }
68 
69 
70 template <typename T>
central_position(const numpy::array_base<T> & array)71 numpy::position central_position(const numpy::array_base<T>& array) {
72     numpy::position res(array.raw_dims(), array.ndims());
73     for (numpy::index_type i = 0, nd = array.ndims(); i != nd; ++i) res.position_[i] /= 2;
74     return res;
75 }
76 
77 template <typename T>
neighbours(const numpy::aligned_array<T> & Bc,bool include_centre=false)78 std::vector<numpy::position> neighbours(const numpy::aligned_array<T>& Bc, bool include_centre = false) {
79     numpy::position centre = central_position(Bc);
80     const unsigned N = Bc.size();
81     typename numpy::aligned_array<T>::const_iterator startc = Bc.begin();
82     std::vector<numpy::position> res;
83     for (unsigned i = 0; i != N; ++i, ++startc) {
84         if (!*startc) continue;
85         if (startc.position() != centre || include_centre) {
86             res.push_back(startc.position() - centre);
87         }
88     }
89     return res;
90 }
91 
92 
93 template <typename T>
neighbours_delta(const numpy::aligned_array<T> & Bc,bool include_centre=false)94 std::vector<numpy::position> neighbours_delta(const numpy::aligned_array<T>& Bc, bool include_centre = false) {
95     std::vector<numpy::position> rs = neighbours(Bc, include_centre);
96     numpy::position accumulated = rs[0];
97     for (unsigned i = 1; i < rs.size(); ++i) {
98         rs[i] -= accumulated;
99         accumulated += rs[i];
100     }
101     return rs;
102 }
103 
104 template<typename T>
margin_of(const numpy::position & position,const numpy::array_base<T> & ref)105 numpy::index_type margin_of(const numpy::position& position, const numpy::array_base<T>& ref) {
106     numpy::index_type margin = std::numeric_limits<numpy::index_type>::max();
107     const numpy::index_type nd = ref.ndims();
108     for (numpy::index_type d = 0; d != nd; ++d) {
109         if (position[d] < margin) margin = position[d];
110         const numpy::index_type rmargin = ref.dim(d) - position[d] - 1;
111         if (rmargin < margin) margin = rmargin;
112    }
113    return margin;
114 }
115 
116 template<typename T>
erode_sub(T a,T b)117 T erode_sub(T a, T b) {
118     if (b == std::numeric_limits<T>::min()) return std::numeric_limits<T>::max();
119     if (!std::numeric_limits<T>::is_signed && (b > a)) return T(0);
120     const T r = a - b;
121     if ( std::numeric_limits<T>::is_signed && (r > a)) return std::numeric_limits<T>::min();
122     return r;
123 }
124 
125 template<>
erode_sub(bool a,bool b)126 bool erode_sub<bool>(bool a, bool b) {
127     return a && b;
128 }
129 
is_bool(T)130 template<typename T> bool is_bool(T) { return false; }
is_bool(bool)131 template<> bool is_bool<bool>(bool) { return true; }
132 
133 template<typename T>
erode(numpy::aligned_array<T> res,const numpy::aligned_array<T> array,const numpy::aligned_array<T> Bc)134 void erode(numpy::aligned_array<T> res, const numpy::aligned_array<T> array, const numpy::aligned_array<T> Bc) {
135     gil_release nogil;
136     const numpy::index_type N = res.size();
137     typename numpy::aligned_array<T>::const_iterator iter = array.begin();
138     filter_iterator<T> filter(array.raw_array(), Bc.raw_array(), ExtendNearest, is_bool(T()));
139     const numpy::index_type N2 = filter.size();
140     T* rpos = res.data();
141     if (!N2) return;
142 
143     for (numpy::index_type i = 0; i != N; ++i, ++rpos, filter.iterate_both(iter)) {
144         T value = std::numeric_limits<T>::max();
145         for (numpy::index_type j = 0; j != N2; ++j) {
146             T arr_val = T();
147             filter.retrieve(iter, j, arr_val);
148             value = std::min<T>(value, erode_sub(arr_val, filter[j]));
149             if (value == std::numeric_limits<T>::min()) break;
150         }
151         *rpos = value;
152     }
153 }
154 
fast_binary_dilate_erode_2d(numpy::aligned_array<bool> res,const numpy::aligned_array<bool> array,numpy::aligned_array<bool> Bc,const bool is_erosion)155 void fast_binary_dilate_erode_2d(numpy::aligned_array<bool> res, const numpy::aligned_array<bool> array, numpy::aligned_array<bool> Bc, const bool is_erosion) {
156     assert(res.is_carray());
157     assert(array.is_carray());
158     assert(res.ndim() == 2);
159 
160     const numpy::index_type Ny = array.dim(0);
161     const numpy::index_type Nx = array.dim(1);
162     const numpy::index_type N = Ny * Nx;
163 
164     const numpy::index_type By = Bc.dim(0);
165     const numpy::index_type Bx = Bc.dim(1);
166 
167     const numpy::index_type Cy = By/2;
168     const numpy::index_type Cx = Bx/2;
169 
170     std::vector<numpy::index_type> positions;
171     for (numpy::index_type y = 0; y != By; ++y) {
172         for (numpy::index_type x = 0; x != Bx; ++x) {
173             if (!Bc.at(y,x)) continue;
174             const numpy::index_type dy = y-Cy;
175             const numpy::index_type dx = x-Cx;
176             if (t_abs(dy) >= Ny || t_abs(dx) >= Nx) continue;
177             if (dy || dx) {
178                 positions.push_back(is_erosion ? dy: -dy);
179                 positions.push_back(is_erosion ? dx: -dx);
180             }
181         }
182     }
183     const numpy::index_type N2 = positions.size()/2;
184     if (Bc.at(Cy,Cx)) std::copy(array.data(), array.data() + N, res.data());
185     else std::fill_n(res.data(), N, is_erosion);
186     if (positions.empty()) return;
187 
188     for (numpy::index_type y = 0; y != Ny; ++y) {
189         bool* const orow = res.data(y);
190         for (numpy::index_type j = 0; j != N2; ++j) {
191             numpy::index_type dy = positions[2*j];
192             numpy::index_type dx = positions[2*j + 1];
193             assert(dx || dy);
194             if ((y + dy) < 0) dy = -y;
195             if ((y + dy) >= Ny) {
196                 dy = -y+(Ny-1);
197             }
198             bool* out = orow;
199             const bool* in = array.data(y + dy);
200             numpy::index_type n = Nx - t_abs(dx);
201             if (dx > 0) {
202                 for (numpy::index_type i = 0; i != (dx-1); ++i) {
203                     if (is_erosion) {
204                         out[Nx-i-1] &= in[Nx-1];
205                     } else {
206                         out[Nx-i-1] |= in[Nx-1];
207                     }
208                 }
209                 in += dx;
210             } else if (dx < 0) {
211                 for (numpy::index_type i = 0; i != -dx; ++i) {
212                     if (is_erosion) {
213                         out[i] &= in[0];
214                     } else {
215                         out[i] |= in[0];
216                     }
217                 }
218                 out += -dx;
219             }
220             if (is_erosion) {
221                 for (numpy::index_type i = 0; i != n; ++i) *out++ &= *in++;
222             } else {
223                 for (numpy::index_type i = 0; i != n; ++i) *out++ |= *in++;
224             }
225         }
226     }
227 }
228 
229 
230 
py_erode(PyObject * self,PyObject * args)231 PyObject* py_erode(PyObject* self, PyObject* args) {
232     PyArrayObject* array;
233     PyArrayObject* Bc;
234     PyArrayObject* output;
235     if (!PyArg_ParseTuple(args, "OOO", &array, &Bc, &output)) return NULL;
236     if (!numpy::are_arrays(array, Bc, output) || !numpy::same_shape(array, output) ||
237         !numpy::equiv_typenums(array, Bc, output) ||
238         PyArray_NDIM(array) != PyArray_NDIM(Bc)
239     ) {
240         PyErr_SetString(PyExc_RuntimeError, TypeErrorMsg);
241         return NULL;
242     }
243     holdref r_o(output);
244 
245     if (numpy::check_type<bool>(array) && PyArray_NDIM(array) == 2 && PyArray_ISCARRAY(array)) {
246         fast_binary_dilate_erode_2d(numpy::aligned_array<bool>(output), numpy::aligned_array<bool>(array), numpy::aligned_array<bool>(Bc), true);
247     } else {
248 #define HANDLE(type) \
249     erode<type>(numpy::aligned_array<type>(output), numpy::aligned_array<type>(array), numpy::aligned_array<type>(Bc));
250     SAFE_SWITCH_ON_INTEGER_TYPES_OF(array);
251 #undef HANDLE
252     }
253 
254     Py_XINCREF(output);
255     return PyArray_Return(output);
256 }
257 
258 template<typename T>
locmin_max(numpy::aligned_array<bool> res,const numpy::aligned_array<T> array,const numpy::aligned_array<T> Bc,bool is_min)259 void locmin_max(numpy::aligned_array<bool> res, const numpy::aligned_array<T> array, const numpy::aligned_array<T> Bc, bool is_min) {
260     gil_release nogil;
261     const numpy::index_type N = res.size();
262     typename numpy::aligned_array<T>::const_iterator iter = array.begin();
263     filter_iterator<T> filter(res.raw_array(), Bc.raw_array(), ExtendNearest, true);
264     const numpy::index_type N2 = filter.size();
265     bool* rpos = res.data();
266 
267     for (numpy::index_type i = 0; i != N; ++i, ++rpos, filter.iterate_both(iter)) {
268         T cur = *iter;
269         for (numpy::index_type j = 0; j != N2; ++j) {
270             T arr_val = T();
271             filter.retrieve(iter, j, arr_val);
272             if (( is_min && (arr_val < cur)) ||
273                 (!is_min && (arr_val > cur))) {
274                     goto skip_to_next;
275                 }
276         }
277         *rpos = true;
278         skip_to_next:
279             ;
280     }
281 }
282 
py_locminmax(PyObject * self,PyObject * args)283 PyObject* py_locminmax(PyObject* self, PyObject* args) {
284     PyArrayObject* array;
285     PyArrayObject* Bc;
286     PyArrayObject* output;
287     int is_min;
288     if (!PyArg_ParseTuple(args, "OOOi", &array, &Bc, &output, &is_min)) return NULL;
289     if (!numpy::are_arrays(array, Bc, output) ||
290         !numpy::same_shape(array, output) ||
291         !numpy::equiv_typenums(array, Bc) ||
292         !numpy::check_type<bool>(output) ||
293         PyArray_NDIM(array) != PyArray_NDIM(Bc) ||
294         !PyArray_ISCARRAY(output)
295     ) {
296         PyErr_SetString(PyExc_RuntimeError, TypeErrorMsg);
297         return NULL;
298     }
299     holdref r_o(output);
300     PyArray_FILLWBYTE(output, 0);
301 
302 #define HANDLE(type) \
303     locmin_max<type>(numpy::aligned_array<bool>(output), numpy::aligned_array<type>(array), numpy::aligned_array<type>(Bc), bool(is_min));
304     SAFE_SWITCH_ON_TYPES_OF(array);
305 #undef HANDLE
306 
307     Py_XINCREF(output);
308     return PyArray_Return(output);
309 }
310 
311 template <typename T>
remove_fake_regmin_max(numpy::aligned_array<bool> regmin,const numpy::aligned_array<T> f,const numpy::aligned_array<T> Bc,bool is_min)312 void remove_fake_regmin_max(numpy::aligned_array<bool> regmin, const numpy::aligned_array<T> f, const numpy::aligned_array<T> Bc, bool is_min) {
313     const numpy::index_type N = f.size();
314     numpy::aligned_array<bool>::iterator riter = regmin.begin();
315     const std::vector<numpy::position> Bc_neighbours = neighbours(Bc);
316     typedef std::vector<numpy::position>::const_iterator Bc_iter;
317     const numpy::index_type N2 = Bc_neighbours.size();
318 
319     for (numpy::index_type i = 0; i != N; ++i, ++riter) {
320         if (!*riter) continue;
321         const numpy::position pos = riter.position();
322         const T val = f.at(pos);
323         for (numpy::index_type j = 0; j != N2; ++j) {
324             numpy::position npos = pos + Bc_neighbours[j];
325             if (f.validposition(npos) &&
326                     !regmin.at(npos) &&
327                             (   (is_min && f.at(npos) <= val) ||
328                                 (!is_min && f.at(npos) >= val)
329                             )) {
330                 numpy::position_stack stack(f.ndims());
331                 assert(regmin.at(pos));
332                 regmin.at(pos) = false;
333                 stack.push(pos);
334                 while (!stack.empty()) {
335                     numpy::position p = stack.top_pop();
336                     for (Bc_iter first = Bc_neighbours.begin(), past = Bc_neighbours.end();
337                                 first != past;
338                                 ++first) {
339                         numpy::position npos = p + *first;
340                         if (regmin.validposition(npos) && regmin.at(npos)) {
341                             regmin.at(npos) = false;
342                             assert(!regmin.at(npos));
343                             stack.push(npos);
344                         }
345                     }
346                 }
347                 // we are done with this position
348                 break;
349             }
350         }
351     }
352 }
353 
py_regminmax(PyObject * self,PyObject * args)354 PyObject* py_regminmax(PyObject* self, PyObject* args) {
355     PyArrayObject* array;
356     PyArrayObject* Bc;
357     PyArrayObject* output;
358     int is_min;
359     if (!PyArg_ParseTuple(args, "OOOi", &array, &Bc, &output, &is_min)) return NULL;
360     if (!numpy::are_arrays(array, Bc, output) || !numpy::same_shape(array, output) ||
361         !PyArray_EquivTypenums(PyArray_TYPE(array), PyArray_TYPE(Bc)) ||
362         !PyArray_EquivTypenums(NPY_BOOL, PyArray_TYPE(output)) ||
363         PyArray_NDIM(array) != PyArray_NDIM(Bc) ||
364         !PyArray_ISCARRAY(output)
365     ) {
366         PyErr_SetString(PyExc_RuntimeError, TypeErrorMsg);
367         return NULL;
368     }
369     holdref r_o(output);
370     PyArray_FILLWBYTE(output, 0);
371 
372 #define HANDLE(type) \
373     locmin_max<type>(numpy::aligned_array<bool>(output), numpy::aligned_array<type>(array), numpy::aligned_array<type>(Bc), bool(is_min)); \
374     remove_fake_regmin_max<type>(numpy::aligned_array<bool>(output), numpy::aligned_array<type>(array), numpy::aligned_array<type>(Bc), bool(is_min));
375 
376     SAFE_SWITCH_ON_TYPES_OF(array);
377 #undef HANDLE
378 
379     Py_XINCREF(output);
380     return PyArray_Return(output);
381 }
382 
383 
384 
385 template <typename T>
dilate_add(T a,T b)386 T dilate_add(T a, T b) {
387     if (a == std::numeric_limits<T>::min()) return a;
388     if (b == std::numeric_limits<T>::min()) return b;
389     const T r = a + b;
390     // if overflow, saturate
391     if (r < std::max<T>(a,b)) return std::numeric_limits<T>::max();
392     return r;
393 }
394 
395 template<>
dilate_add(bool a,bool b)396 bool dilate_add(bool a, bool b) {
397     return a && b;
398 }
399 
400 template<typename T>
dilate(numpy::aligned_array<T> res,const numpy::array<T> array,const numpy::aligned_array<T> Bc)401 void dilate(numpy::aligned_array<T> res, const numpy::array<T> array, const numpy::aligned_array<T> Bc) {
402     gil_release nogil;
403     const numpy::index_type N = res.size();
404     typename numpy::array<T>::const_iterator iter = array.begin();
405     filter_iterator<T> filter(res.raw_array(), Bc.raw_array(), ExtendNearest, is_bool(T()));
406     const numpy::index_type N2 = filter.size();
407     // T* is a fine iterator type.
408     T* rpos = res.data();
409     std::fill(rpos, rpos + res.size(), std::numeric_limits<T>::min());
410     if (!N2) return;
411 
412     for (numpy::index_type i = 0; i != N; ++i, ++rpos, filter.iterate_both(iter)) {
413         const T value = *iter;
414         if (value == std::numeric_limits<T>::min()) continue;
415         for (numpy::index_type j = 0; j != N2; ++j) {
416             const T nval = dilate_add(value, filter[j]);
417             T arr_val = T();
418             filter.retrieve(rpos, j, arr_val);
419             if (nval > arr_val) filter.set(rpos, j, nval);
420         }
421     }
422 }
423 
424 
py_dilate(PyObject * self,PyObject * args)425 PyObject* py_dilate(PyObject* self, PyObject* args) {
426     PyArrayObject* array;
427     PyArrayObject* Bc;
428     PyArrayObject* output;
429     if (!PyArg_ParseTuple(args,"OOO", &array, &Bc, &output)) return NULL;
430     if (!numpy::are_arrays(array, Bc, output) ||
431         !numpy::same_shape(array, output) ||
432         !numpy::equiv_typenums(array, Bc, output) ||
433         PyArray_NDIM(array) != PyArray_NDIM(Bc)
434     ) {
435         PyErr_SetString(PyExc_RuntimeError, TypeErrorMsg);
436         return NULL;
437     }
438     holdref r_o(output);
439     if (numpy::check_type<bool>(array) && PyArray_NDIM(array) == 2 && PyArray_ISCARRAY(array)) {
440         fast_binary_dilate_erode_2d(numpy::aligned_array<bool>(output), numpy::aligned_array<bool>(array), numpy::aligned_array<bool>(Bc), false);
441     } else {
442 #define HANDLE(type) \
443         dilate<type>(numpy::aligned_array<type>(output),numpy::array<type>(array),numpy::aligned_array<type>(Bc));
444         SAFE_SWITCH_ON_INTEGER_TYPES_OF(array);
445 #undef HANDLE
446     }
447 
448     Py_XINCREF(output);
449     return PyArray_Return(output);
450 }
451 
py_disk_2d(PyObject * self,PyObject * args)452 PyObject* py_disk_2d(PyObject* self, PyObject* args) {
453     PyArrayObject* array;
454     int radius;
455     if (!PyArg_ParseTuple(args,"Oi", &array, &radius)) return NULL;
456     if (!numpy::are_arrays(array) ||
457         PyArray_NDIM(array) != 2 ||
458         !PyArray_ISCARRAY(array) ||
459         !numpy::check_type<bool>(array) ||
460         radius < 0
461     ) {
462         PyErr_SetString(PyExc_RuntimeError, TypeErrorMsg);
463         return NULL;
464     }
465     Py_XINCREF(array);
466     bool* iter = numpy::ndarray_cast<bool*>(array);
467     const int radius2 = radius * radius;
468     const numpy::index_type N0 = PyArray_DIM(array, 0);
469     const numpy::index_type N1 = PyArray_DIM(array, 1);
470     const numpy::index_type c0 = N0/2;
471     const numpy::index_type c1 = N1/2;
472     for (numpy::index_type x0 = 0; x0 != N0; ++x0) {
473         for (numpy::index_type x1 = 0; x1 != N1; ++x1, ++iter) {
474             if ((x0-c0)*(x0-c0) + (x1-c1)*(x1-c1) < radius2) {
475                 *iter = true;
476             }
477         }
478     }
479     return PyArray_Return(array);
480 }
481 
close_holes(const numpy::aligned_array<bool> ref,numpy::aligned_array<bool> f,const numpy::aligned_array<bool> Bc)482 void close_holes(const numpy::aligned_array<bool> ref, numpy::aligned_array<bool> f, const numpy::aligned_array<bool> Bc) {
483     std::fill_n(f.data(), f.size(), false);
484 
485     numpy::position_stack stack(ref.ndim());
486     const numpy::index_type N = ref.size();
487     const std::vector<numpy::position> Bc_neighbours = neighbours(Bc);
488     const numpy::index_type N2 = Bc_neighbours.size();
489     for (numpy::index_type d = 0; d != ref.ndims(); ++d) {
490         if (ref.dim(d) == 0) continue;
491         numpy::position pos;
492         pos.nd_ = ref.ndims();
493         for (numpy::index_type di = 0; di != ref.ndims(); ++di) pos.position_[di] = 0;
494 
495         for (numpy::index_type i = 0; i != N/ref.dim(d); ++i) {
496             pos.position_[d] = 0;
497             if (!ref.at(pos) && !f.at(pos)) {
498                 f.at(pos) = true;
499                 stack.push(pos);
500             }
501             pos.position_[d] = ref.dim(d) - 1;
502             if (!ref.at(pos) && !f.at(pos)) {
503                 f.at(pos) = true;
504                 stack.push(pos);
505             }
506 
507             for (numpy::index_type j = 0; j != ref.ndims() - 1; ++j) {
508                 if (j == d) ++j;
509                 if (pos.position_[j] < numpy::index_type(ref.dim(j))) {
510                     ++pos.position_[j];
511                     break;
512                 }
513                 pos.position_[j] = 0;
514             }
515         }
516     }
517     while (!stack.empty()) {
518         numpy::position pos = stack.top_pop();
519         std::vector<numpy::position>::const_iterator startc = Bc_neighbours.begin();
520         for (numpy::index_type j = 0; j != N2; ++j, ++startc) {
521             numpy::position npos = pos + *startc;
522             if (ref.validposition(npos) && !ref.at(npos) && !f.at(npos)) {
523                 f.at(npos) = true;
524                 stack.push(npos);
525             }
526         }
527     }
528     for (bool* start = f.data(), *past = f.data() + f.size(); start != past; ++ start) {
529         *start = !*start;
530     }
531 }
532 
py_close_holes(PyObject * self,PyObject * args)533 PyObject* py_close_holes(PyObject* self, PyObject* args) {
534     PyArrayObject* ref;
535     PyArrayObject* Bc;
536     if (!PyArg_ParseTuple(args,"OO", &ref, &Bc)) return NULL;
537     if (!numpy::are_arrays(ref, Bc) ||
538         !numpy::equiv_typenums(ref, Bc) ||
539         !numpy::check_type<bool>(ref)
540         ) {
541         PyErr_SetString(PyExc_RuntimeError,TypeErrorMsg);
542         return NULL;
543     }
544     PyArrayObject* res_a = (PyArrayObject*)PyArray_SimpleNew(PyArray_NDIM(ref),
545                                             PyArray_DIMS(ref),
546                                             PyArray_TYPE(ref));
547     if (!res_a) return NULL;
548 
549     // We own the only reference.
550     // If an exception happens, we want r_a to delete the object.
551     // So we do call it with incref=false:
552     holdref r_a(res_a, false);
553     try {
554         close_holes(numpy::aligned_array<bool>(ref), numpy::aligned_array<bool>(res_a), numpy::aligned_array<bool>(Bc));
555     }
556     CATCH_PYTHON_EXCEPTIONS
557 
558     Py_INCREF(res_a);
559     return PyArray_Return(res_a);
560 }
561 
562 template <typename CostType>
563 struct MarkerInfo {
564     CostType cost;
565     npy_intp idx;
566     npy_intp position;
567     npy_intp margin;
MarkerInfo__anonc362df840111::MarkerInfo568     MarkerInfo(CostType cost, npy_intp idx, npy_intp position, npy_intp margin)
569         :cost(cost)
570         ,idx(idx)
571         ,position(position)
572         ,margin(margin) {
573         }
operator <__anonc362df840111::MarkerInfo574     bool operator < (const MarkerInfo& other) const {
575         // We want the smaller argument to compare higher, so we reverse the order here:
576         if (cost == other.cost) return idx > other.idx;
577         return cost > other.cost;
578     }
579 };
580 
581 struct NeighbourElem {
582     // This represents a neighbour
583     //  - delta: how to get from the current element to the next by pointer manipulation
584     //  - step: how large the distance to the centre is
585     //  - delta_position: the difference as a numpy::position
NeighbourElem__anonc362df840111::NeighbourElem586     NeighbourElem(npy_intp delta, npy_intp step, const numpy::position& delta_position)
587         :delta(delta)
588         ,step(step)
589         ,delta_position(delta_position)
590         { }
591     npy_intp delta;
592     npy_intp step;
593     numpy::position delta_position;
594 };
595 
596 template<typename BaseType>
cwatershed(numpy::aligned_array<npy_int64> res,numpy::aligned_array<bool> * lines,const numpy::aligned_array<BaseType> array,const numpy::aligned_array<npy_int64> markers,const numpy::aligned_array<BaseType> Bc)597 void cwatershed(numpy::aligned_array<npy_int64> res,
598                         numpy::aligned_array<bool>* lines,
599                         const numpy::aligned_array<BaseType> array,
600                         const numpy::aligned_array<npy_int64> markers,
601                         const numpy::aligned_array<BaseType> Bc) {
602     gil_release nogil;
603     const npy_intp N = res.size();
604     const npy_intp N2 = Bc.size();
605     assert(res.is_carray());
606     npy_int64* rdata = res.data();
607     std::vector<NeighbourElem> neighbours;
608     const numpy::position centre = central_position(Bc);
609     typename numpy::aligned_array<BaseType>::const_iterator Bi = Bc.begin();
610     for (npy_intp j = 0; j != N2; ++j, ++Bi) {
611         if (*Bi) {
612             numpy::position npos = Bi.position() - centre;
613             npy_intp margin = 0;
614             for (numpy::index_type d = 0; d != Bc.ndims(); ++d) {
615                 margin = std::max<npy_intp>(t_abs(npy_intp(npos[d])), margin);
616             }
617             npy_intp delta = markers.pos_to_flat(npos);
618             if (!delta) continue;
619             neighbours.push_back(NeighbourElem(delta, margin, npos));
620         }
621     }
622     npy_intp idx = 0;
623 
624     enum { white, grey, black };
625     std::vector<unsigned char> status(array.size(), static_cast<unsigned char>(white));
626 
627     std::priority_queue<MarkerInfo<BaseType> > hqueue;
628 
629     typename numpy::aligned_array<npy_int64>::const_iterator miter = markers.begin();
630     for (numpy::index_type i = 0; i != N; ++i, ++miter) {
631         if (*miter) {
632             assert(markers.validposition(miter.position()));
633             const numpy::position mpos = miter.position();
634             const npy_intp margin = margin_of(mpos, markers);
635 
636             hqueue.push(MarkerInfo<BaseType>(array.at(mpos), idx++, markers.pos_to_flat(mpos), margin));
637             res.at(mpos) = *miter;
638             status[markers.pos_to_flat(mpos)] = grey;
639         }
640     }
641 
642     while (!hqueue.empty()) {
643         const MarkerInfo<BaseType> next = hqueue.top();
644         hqueue.pop();
645         status[next.position] = black;
646         npy_intp margin = next.margin;
647         for (typename std::vector<NeighbourElem>::const_iterator neighbour = neighbours.begin(),
648                             past = neighbours.end();
649                     neighbour != past;
650                     ++neighbour) {
651             const numpy::index_type npos = next.position + neighbour->delta;
652             numpy::index_type nmargin = margin - neighbour->step;
653             if (nmargin < 0) {
654                 // nmargin is a lower bound on the margin, so we must recompute the actual value
655                 numpy::position pos = markers.flat_to_pos(next.position);
656                 assert(markers.validposition(pos));
657                 numpy::position long_pos = pos + neighbour->delta_position;
658                 nmargin = margin_of(long_pos, markers);
659                 if (nmargin < 0) {
660                     // We are outside the image.
661                     continue;
662                 }
663                 // we are good with the recomputed margin
664                 assert(markers.validposition(long_pos));
665                 // Update lower bound
666                 if ((nmargin - neighbour->step) > margin) margin = nmargin - neighbour->step;
667             }
668             assert(npos < npy_intp(status.size()));
669             switch (status[npos]) {
670                 case white: {
671                     const BaseType ncost = array.at_flat(npos);
672                     rdata[npos] = rdata[next.position];
673                     hqueue.push(MarkerInfo<BaseType>(ncost, idx++, npos, nmargin));
674                     status[npos] = grey;
675                     break;
676                 }
677                 case grey: {
678                     if (lines && rdata[next.position] != rdata[npos]) {
679                         lines->at_flat(npos) = true;
680                     }
681                     break;
682                 }
683             }
684         }
685     }
686 }
687 
py_cwatershed(PyObject * self,PyObject * args)688 PyObject* py_cwatershed(PyObject* self, PyObject* args) {
689     PyArrayObject* array;
690     PyArrayObject* markers;
691     PyArrayObject* Bc;
692     int return_lines;
693     if (!PyArg_ParseTuple(args,"OOOi", &array, &markers, &Bc,&return_lines)) {
694         return NULL;
695     }
696     if (!numpy::are_arrays(array, markers, Bc) ||
697         !numpy::check_type<npy_int64>(markers)) {
698         PyErr_SetString(PyExc_RuntimeError, "mahotas._cwatershed: markers should be an int32 array.");
699         return NULL;
700     }
701     PyArrayObject* res_a = (PyArrayObject*)PyArray_SimpleNew(
702                                                     PyArray_NDIM(array),
703                                                     PyArray_DIMS(array),
704                                                     NPY_INT64);
705     if (!res_a) return NULL;
706     PyArrayObject* lines =  0;
707     numpy::aligned_array<bool>* lines_a = 0;
708     if (return_lines) {
709         lines = (PyArrayObject*)PyArray_SimpleNew(PyArray_NDIM(array), PyArray_DIMS(array), NPY_BOOL);
710         if (!lines) return NULL;
711         lines_a = new numpy::aligned_array<bool>(lines);
712     }
713 #define HANDLE(type) \
714     cwatershed<type>(numpy::aligned_array<npy_int64>(res_a),lines_a,numpy::aligned_array<type>(array),numpy::aligned_array<npy_int64>(markers),numpy::aligned_array<type>(Bc));
715     SAFE_SWITCH_ON_TYPES_OF(array);
716 #undef HANDLE
717     if (return_lines) {
718         delete lines_a;
719         PyObject* ret_val = PyTuple_New(2);
720         PyTuple_SetItem(ret_val,0,(PyObject*)res_a);
721         PyTuple_SetItem(ret_val,1,(PyObject*)lines);
722         return ret_val;
723     }
724     return PyArray_Return(res_a);
725 }
726 
727 
compute_euc2_dist(const numpy::position & a,const numpy::position b)728 double compute_euc2_dist(const numpy:: position& a, const numpy::position b) {
729     assert(a.ndim() == b.ndim());
730     double r = 0.;
731     const numpy::index_type n = a.ndim();
732     for (numpy::index_type i = 0; i != n; ++i) {
733         r += (a[i]-b[i])*(a[i]-b[i]);
734     }
735     return r;
736 }
737 
738 
739 // Arguably the function distance_multi should have been in file distance.cpp,
740 // but it uses several functions from this module like neighbours().
741 template<typename BaseType>
distance_multi(numpy::aligned_array<BaseType> res,const numpy::aligned_array<bool> array,const numpy::aligned_array<bool> Bc)742 void distance_multi(numpy::aligned_array<BaseType> res,
743                         const numpy::aligned_array<bool> array,
744                         const numpy::aligned_array<bool> Bc) {
745     gil_release nogil;
746     const numpy::index_type N = res.size();
747     const std::vector<numpy::position> Bcs = neighbours_delta(Bc);
748     const numpy::index_type N2 = Bcs.size();
749 
750     typename numpy::aligned_array<bool>::const_iterator aiter = array.begin();
751     typename numpy::aligned_array<BaseType>::iterator riter = res.begin();
752 
753     numpy::position_queue cur_q(res.ndim());
754     numpy::position_queue orig_q(res.ndim());
755     std::queue<double> dist_q;
756     for (numpy::index_type i = 0; i != N; ++i, ++riter, ++aiter) {
757         if (!*aiter) {
758             *riter = 0;
759             const numpy::position p = aiter.position();
760             numpy::position next = p;
761             for (numpy::index_type j = 0; j != N2; ++j) {
762                 next += Bcs[j];
763                 if (array.validposition(next) && array.at(next)) {
764                     const double dist = compute_euc2_dist(next, p);
765                     BaseType* rpos = res.data(next);
766                     if (*rpos > dist) {
767                         *rpos = dist;
768 
769                         cur_q.push(next);
770                         orig_q.push(p);
771                         dist_q.push(dist);
772                     }
773                 }
774             }
775         }
776     }
777 
778     while (!dist_q.empty()) {
779         numpy::position next = cur_q.top_pop();
780         const numpy::position orig = orig_q.top_pop();
781         const BaseType dist = dist_q.front();
782         dist_q.pop();
783 
784         assert(dist == compute_euc2_dist(next, orig));
785 
786         if (res.at(next) < dist) continue;
787         for (numpy::index_type j = 0; j != N2; ++j) {
788             next += Bcs[j];
789             if (array.validposition(next)) {
790                 const double next_dist = compute_euc2_dist(next, orig);
791                 BaseType* rpos = res.data(next);
792                 if (*rpos > next_dist) {
793                     *rpos = next_dist;
794                     cur_q.push(next);
795                     orig_q.push(orig);
796                     dist_q.push(next_dist);
797                 }
798             }
799         }
800     }
801 }
802 
py_distance_multi(PyObject * self,PyObject * args)803 PyObject* py_distance_multi(PyObject* self, PyObject* args) {
804     PyArrayObject* array;
805     PyArrayObject* res;
806     PyArrayObject* Bc;
807     if (!PyArg_ParseTuple(args,"OOO", &res, &array, &Bc)) {
808         return NULL;
809     }
810     if (!numpy::are_arrays(array, res, Bc) ||
811         !numpy::check_type<bool>(array) ||
812         !numpy::check_type<bool>(Bc) ||
813         !numpy::same_shape(array, res)
814         ) {
815         PyErr_SetString(PyExc_RuntimeError, "mahotas._distance_multi: res and input array should have same shape. input & Bc arrays maust be boolean arrays.");
816         return NULL;
817     }
818 #define HANDLE(type) \
819     distance_multi<type>(numpy::aligned_array<type>(res), \
820                         numpy::aligned_array<bool>(array),\
821                         numpy::aligned_array<bool>(Bc));
822 
823     SAFE_SWITCH_ON_TYPES_OF(res);
824 #undef HANDLE
825 
826     Py_RETURN_NONE;
827 }
828 
829 struct HitMissNeighbour {
HitMissNeighbour__anonc362df840111::HitMissNeighbour830     HitMissNeighbour(numpy::index_type delta, int value)
831         :delta(delta)
832         ,value(value)
833         { }
834     numpy::index_type delta;
835     int value;
836 };
837 
838 template <typename T>
hitmiss(numpy::aligned_array<T> res,const numpy::aligned_array<T> & input,const numpy::aligned_array<T> & Bc)839 void hitmiss(numpy::aligned_array<T> res, const numpy::aligned_array<T>& input, const numpy::aligned_array<T>& Bc) {
840     gil_release nogil;
841     typedef typename numpy::aligned_array<T>::const_iterator const_iterator;
842     const numpy::index_type N = input.size();
843     const numpy::index_type N2 = Bc.size();
844     const numpy::position centre = central_position(Bc);
845     numpy::index_type Bc_margin = 0;
846     for (numpy::index_type d = 0; d != Bc.ndims(); ++d) {
847         numpy::index_type cmargin = Bc.dim(d)/2;
848         if (cmargin > Bc_margin) Bc_margin = cmargin;
849     }
850 
851     std::vector<HitMissNeighbour> neighbours;
852     const_iterator Bi = Bc.begin();
853     for (numpy::index_type j = 0; j != N2; ++j, ++Bi) {
854         if (*Bi != 2) {
855             numpy::position npos = Bi.position() - centre;
856             numpy::index_type delta = input.pos_to_flat(npos);
857             neighbours.push_back(HitMissNeighbour(delta, *Bi));
858         }
859     }
860 
861     // This is a micro-optimisation for templates with structure
862     // It makes it more likely that matching will fail earlier
863     // in uniform regions than otherwise would be the case.
864     std::random_shuffle(neighbours.begin(), neighbours.end());
865     numpy::index_type slack = 0;
866     for (npy_intp i = 0; i != N; ++i) {
867         while (!slack) {
868             numpy::position cur = input.flat_to_pos(i);
869             bool moved = false;
870             for (numpy::index_type d = 0; d != input.ndims(); ++d) {
871                 numpy::index_type margin = std::min<numpy::index_type>(cur[d], input.dim(d) - cur[d] - 1);
872                 if (margin < Bc.dim(d)/2) {
873                     numpy::index_type size = 1;
874                     for (numpy::index_type dd = d+1; dd < input.ndims(); ++dd) size *= input.dim(dd);
875                     for (numpy::index_type j = 0; j != size; ++j) {
876                         res.at_flat(i++) = 0;
877                         if (i == N) return;
878                     }
879                     moved = true;
880                     break;
881                 }
882             }
883             if (!moved) slack = input.dim(input.ndims() - 1) - Bc.dim(input.ndims() - 1) + 1;
884         }
885         --slack;
886         T value = 1;
887         for (std::vector<HitMissNeighbour>::const_iterator neighbour = neighbours.begin(), past = neighbours.end();
888             neighbour != past;
889             ++neighbour) {
890             if (input.at_flat(i + neighbour->delta) != static_cast<T>(neighbour->value)) {
891                 value = 0;
892                 break;
893             }
894         }
895         res.at_flat(i) = value;
896     }
897 }
898 
py_hitmiss(PyObject * self,PyObject * args)899 PyObject* py_hitmiss(PyObject* self, PyObject* args) {
900     PyArrayObject* array;
901     PyArrayObject* Bc;
902     PyArrayObject* res_a;
903     if (!PyArg_ParseTuple(args, "OOO", &array, &Bc, &res_a)) {
904         PyErr_SetString(PyExc_RuntimeError,TypeErrorMsg);
905         return NULL;
906     }
907     holdref r(res_a);
908 
909 #define HANDLE(type) \
910     hitmiss<type>(numpy::aligned_array<type>(res_a), numpy::aligned_array<type>(array), numpy::aligned_array<type>(Bc));
911     SAFE_SWITCH_ON_INTEGER_TYPES_OF(array);
912 #undef HANDLE
913 
914     Py_INCREF(res_a);
915     return PyArray_Return(res_a);
916 }
917 
py_majority_filter(PyObject * self,PyObject * args)918 PyObject* py_majority_filter(PyObject* self, PyObject* args) {
919     PyArrayObject* array;
920     PyArrayObject* res_a;
921     long long N;
922     if (!PyArg_ParseTuple(args, "OLO", &array, &N, &res_a) ||
923         !PyArray_Check(array) || !PyArray_Check(res_a) ||
924         PyArray_TYPE(array) != NPY_BOOL || PyArray_TYPE(res_a) != NPY_BOOL ||
925         !PyArray_ISCARRAY(res_a)) {
926         PyErr_SetString(PyExc_RuntimeError,TypeErrorMsg);
927         return NULL;
928     }
929     Py_INCREF(res_a);
930     PyArray_FILLWBYTE(res_a, 0);
931     numpy::aligned_array<bool> input(array);
932     numpy::aligned_array<bool> output(res_a);
933     const numpy::index_type rows = input.dim(0);
934     const numpy::index_type cols = input.dim(1);
935     const numpy::index_type T = N*N/2;
936     if (rows < N || cols < N) {
937         return PyArray_Return(res_a);
938     }
939     for (numpy::index_type y = 0; y != rows-N; ++y) {
940         bool* output_iter = output.data() + (y+numpy::index_type(N/2)) * output.stride(0) + numpy::index_type(N/2);
941         for (numpy::index_type x = 0; x != cols-N; ++x) {
942             numpy::index_type count = 0;
943             for (numpy::index_type dy = 0; dy != N; ++dy) {
944                 for (numpy::index_type dx = 0; dx != N; ++dx) {
945                     if (input.at(y+dy,x+dx)) ++count;
946                 }
947             }
948             if (count >= T) {
949                 *output_iter = true;
950             }
951             ++output_iter;
952         }
953     }
954     return PyArray_Return(res_a);
955 }
956 
957 
958 
959 
960 PyMethodDef methods[] = {
961   {"subm",(PyCFunction)py_subm, METH_VARARGS, NULL},
962   {"dilate",(PyCFunction)py_dilate, METH_VARARGS, NULL},
963   {"disk_2d",(PyCFunction)py_disk_2d, METH_VARARGS, NULL},
964   {"erode",(PyCFunction)py_erode, METH_VARARGS, NULL},
965   {"close_holes",(PyCFunction)py_close_holes, METH_VARARGS, NULL},
966   {"cwatershed",(PyCFunction)py_cwatershed, METH_VARARGS, NULL},
967   {"distance_multi",(PyCFunction)py_distance_multi, METH_VARARGS, NULL},
968   {"locmin_max",(PyCFunction)py_locminmax, METH_VARARGS, NULL},
969   {"regmin_max",(PyCFunction)py_regminmax, METH_VARARGS, NULL},
970   {"hitmiss",(PyCFunction)py_hitmiss, METH_VARARGS, NULL},
971   {"majority_filter",(PyCFunction)py_majority_filter, METH_VARARGS, NULL},
972   {NULL, NULL,0,NULL},
973 };
974 
975 } // namespace
976 
977 DECLARE_MODULE(_morph)
978