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