1#cython: cdivision=True
2#cython: boundscheck=False
3#cython: nonecheck=False
4#cython: wraparound=False
5
6import numpy as np
7
8cimport numpy as cnp
9from libc.stdlib cimport malloc, free
10
11cnp.import_array()
12
13cdef inline dtype_t _max(dtype_t a, dtype_t b) nogil:
14    return a if a >= b else b
15
16
17cdef inline dtype_t _min(dtype_t a, dtype_t b) nogil:
18    return a if a <= b else b
19
20
21cdef inline void histogram_increment(Py_ssize_t[::1] histo, double* pop,
22                                     dtype_t value) nogil:
23    histo[value] += 1
24    pop[0] += 1
25
26
27cdef inline void histogram_decrement(Py_ssize_t[::1] histo, double* pop,
28                                     dtype_t value) nogil:
29    histo[value] -= 1
30    pop[0] -= 1
31
32
33cdef inline char is_in_mask(Py_ssize_t rows, Py_ssize_t cols,
34                            Py_ssize_t r, Py_ssize_t c,
35                            char* mask) nogil:
36    """Check whether given coordinate is within image and mask is true."""
37    if r < 0 or r > rows - 1 or c < 0 or c > cols - 1:
38        return 0
39    else:
40        if not mask:
41            return 1
42        elif mask[r * cols + c]:
43            return 1
44        else:
45            return 0
46
47
48cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t[::1], double,
49                            dtype_t, Py_ssize_t, Py_ssize_t, double,
50                            double, Py_ssize_t, Py_ssize_t) nogil,
51                dtype_t[:, ::1] image,
52                char[:, ::1] footprint,
53                char[:, ::1] mask,
54                dtype_t_out[:, :, ::1] out,
55                signed char shift_x, signed char shift_y,
56                double p0, double p1,
57                Py_ssize_t s0, Py_ssize_t s1,
58                Py_ssize_t n_bins) except *:
59    """Compute histogram for each pixel neighborhood, apply kernel function and
60    use kernel function return value for output image.
61    """
62
63    cdef Py_ssize_t rows = image.shape[0]
64    cdef Py_ssize_t cols = image.shape[1]
65    cdef Py_ssize_t srows = footprint.shape[0]
66    cdef Py_ssize_t scols = footprint.shape[1]
67    cdef Py_ssize_t odepth = out.shape[2]
68
69    cdef Py_ssize_t centre_r = <Py_ssize_t>(footprint.shape[0] / 2) + shift_y
70    cdef Py_ssize_t centre_c = <Py_ssize_t>(footprint.shape[1] / 2) + shift_x
71
72    # check that footprint center is inside the element bounding box
73    assert centre_r >= 0, f'centre_r {centre_r} < 0'
74    assert centre_c >= 0, f'centre_c {centre_c} < 0'
75    assert centre_r < srows, f'centre_r {centre_r} >= srows {srows}'
76    assert centre_c < scols, f'centre_c {centre_c} >= scols {scols}'
77
78    cdef Py_ssize_t mid_bin = n_bins / 2
79
80    # define pointers to the data
81    cdef char * mask_data = NULL
82    if mask is not None:
83        mask_data = &mask[0, 0]
84
85    # define local variable types
86    cdef Py_ssize_t r, c, rr, cc, s, value, local_max, i, even_row
87
88    # number of pixels actually inside the neighborhood (double)
89    cdef double pop = 0
90
91    # build attack and release borders by using difference along axis
92    t = np.hstack((footprint, np.zeros((footprint.shape[0], 1))))
93    cdef unsigned char[:, :] t_e = (np.diff(t, axis=1) < 0).view(np.uint8)
94
95    t = np.hstack((np.zeros((footprint.shape[0], 1)), footprint))
96    cdef unsigned char[:, :] t_w = (np.diff(t, axis=1) > 0).view(np.uint8)
97
98    t = np.vstack((footprint, np.zeros((1, footprint.shape[1]))))
99    cdef unsigned char[:, :] t_s = (np.diff(t, axis=0) < 0).view(np.uint8)
100
101    t = np.vstack((np.zeros((1, footprint.shape[1])), footprint))
102    cdef unsigned char[:, :] t_n = (np.diff(t, axis=0) > 0).view(np.uint8)
103
104    # the current local histogram distribution
105    # cdef Py_ssize_t* histo
106    cdef Py_ssize_t [::1] histo = np.zeros(n_bins, dtype=np.intp)
107
108    # these lists contain the relative pixel row and column for each of the 4
109    # attack borders east, west, north and south e.g. se_e_r lists the rows of
110    # the east footprint border
111    cdef Py_ssize_t se_size = srows * scols
112    cdef Py_ssize_t [::1] se_e_r = np.empty(se_size, dtype=np.intp)
113    cdef Py_ssize_t [::1] se_e_c = np.empty(se_size, dtype=np.intp)
114
115    cdef Py_ssize_t [::1] se_w_r = np.empty(se_size, dtype=np.intp)
116    cdef Py_ssize_t [::1] se_w_c = np.empty(se_size, dtype=np.intp)
117
118    cdef Py_ssize_t [::1] se_n_r = np.empty(se_size, dtype=np.intp)
119    cdef Py_ssize_t [::1] se_n_c = np.empty(se_size, dtype=np.intp)
120
121    cdef Py_ssize_t [::1] se_s_r = np.empty(se_size, dtype=np.intp)
122    cdef Py_ssize_t [::1] se_s_c = np.empty(se_size, dtype=np.intp)
123
124    # number of element in each attack border
125    cdef Py_ssize_t num_se_n, num_se_s, num_se_e, num_se_w
126
127    num_se_n = num_se_s = num_se_e = num_se_w = 0
128
129    with nogil:
130        for i in range(n_bins):
131            histo[i] = 0
132
133        for r in range(srows):
134            for c in range(scols):
135                if t_e[r, c]:
136                    se_e_r[num_se_e] = r - centre_r
137                    se_e_c[num_se_e] = c - centre_c
138                    num_se_e += 1
139                if t_w[r, c]:
140                    se_w_r[num_se_w] = r - centre_r
141                    se_w_c[num_se_w] = c - centre_c
142                    num_se_w += 1
143                if t_n[r, c]:
144                    se_n_r[num_se_n] = r - centre_r
145                    se_n_c[num_se_n] = c - centre_c
146                    num_se_n += 1
147                if t_s[r, c]:
148                    se_s_r[num_se_s] = r - centre_r
149                    se_s_c[num_se_s] = c - centre_c
150                    num_se_s += 1
151
152        for r in range(srows):
153            for c in range(scols):
154                rr = r - centre_r
155                cc = c - centre_c
156                if footprint[r, c]:
157                    if is_in_mask(rows, cols, rr, cc, mask_data):
158                        histogram_increment(histo, &pop, image[rr, cc])
159
160        r = 0
161        c = 0
162        kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], n_bins, mid_bin,
163               p0, p1, s0, s1)
164
165        # main loop
166        r = 0
167        for even_row in range(0, rows, 2):
168
169            # ---> west to east
170            for c in range(1, cols):
171                for s in range(num_se_e):
172                    rr = r + se_e_r[s]
173                    cc = c + se_e_c[s]
174                    if is_in_mask(rows, cols, rr, cc, mask_data):
175                        histogram_increment(histo, &pop, image[rr, cc])
176
177                for s in range(num_se_w):
178                    rr = r + se_w_r[s]
179                    cc = c + se_w_c[s] - 1
180                    if is_in_mask(rows, cols, rr, cc, mask_data):
181                        histogram_decrement(histo, &pop, image[rr, cc])
182
183                kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], n_bins,
184                       mid_bin, p0, p1, s0, s1)
185
186            r += 1  # pass to the next row
187            if r >= rows:
188                break
189
190            # ---> north to south
191            for s in range(num_se_s):
192                rr = r + se_s_r[s]
193                cc = c + se_s_c[s]
194                if is_in_mask(rows, cols, rr, cc, mask_data):
195                    histogram_increment(histo, &pop, image[rr, cc])
196
197            for s in range(num_se_n):
198                rr = r + se_n_r[s] - 1
199                cc = c + se_n_c[s]
200                if is_in_mask(rows, cols, rr, cc, mask_data):
201                    histogram_decrement(histo, &pop, image[rr, cc])
202
203            kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], n_bins,
204                   mid_bin, p0, p1, s0, s1)
205
206            # ---> east to west
207            for c in range(cols - 2, -1, -1):
208                for s in range(num_se_w):
209                    rr = r + se_w_r[s]
210                    cc = c + se_w_c[s]
211                    if is_in_mask(rows, cols, rr, cc, mask_data):
212                        histogram_increment(histo, &pop, image[rr, cc])
213
214                for s in range(num_se_e):
215                    rr = r + se_e_r[s]
216                    cc = c + se_e_c[s] + 1
217                    if is_in_mask(rows, cols, rr, cc, mask_data):
218                        histogram_decrement(histo, &pop, image[rr, cc])
219
220                kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], n_bins,
221                       mid_bin, p0, p1, s0, s1)
222
223            r += 1  # pass to the next row
224            if r >= rows:
225                break
226
227            # ---> north to south
228            for s in range(num_se_s):
229                rr = r + se_s_r[s]
230                cc = c + se_s_c[s]
231                if is_in_mask(rows, cols, rr, cc, mask_data):
232                    histogram_increment(histo, &pop, image[rr, cc])
233
234            for s in range(num_se_n):
235                rr = r + se_n_r[s] - 1
236                cc = c + se_n_c[s]
237                if is_in_mask(rows, cols, rr, cc, mask_data):
238                    histogram_decrement(histo, &pop, image[rr, cc])
239
240            kernel(&out[r, c, 0], odepth, histo, pop, image[r, c],
241                   n_bins, mid_bin, p0, p1, s0, s1)
242