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