scikit-image

1
#cython: cdivision=True
2
#cython: boundscheck=False
3
#cython: nonecheck=False
4
#cython: wraparound=False
5
import numpy as np
6
cimport numpy as cnp
7
from libc.float cimport DBL_MAX
8
from libc.math cimport atan2, fabs
9

10
from .._shared.fused_numerics cimport np_floats
11

12
cnp.import_array()
13

14

15
def _corner_moravec(np_floats[:, ::1] cimage, Py_ssize_t window_size=1):
16
    """Compute Moravec corner measure response image.
17

18
    This is one of the simplest corner detectors and is comparatively fast but
19
    has several limitations (e.g. not rotation invariant).
20

21
    Parameters
22
    ----------
23
    image : ndarray
24
        Input image.
25
    window_size : int, optional (default 1)
26
        Window size.
27

28
    Returns
29
    -------
30
    response : ndarray
31
        Moravec response image.
32

33
    References
34
    ----------
35
    .. [1] http://kiwi.cs.dal.ca/~dparks/CornerDetection/moravec.htm
36
    .. [2] https://en.wikipedia.org/wiki/Corner_detection
37

38
    Examples
39
    --------
40
    >>> from skimage.feature import corner_moravec
41
    >>> square = np.zeros([7, 7])
42
    >>> square[3, 3] = 1
43
    >>> square.astype(int)
44
    array([[0, 0, 0, 0, 0, 0, 0],
45
           [0, 0, 0, 0, 0, 0, 0],
46
           [0, 0, 0, 0, 0, 0, 0],
47
           [0, 0, 0, 1, 0, 0, 0],
48
           [0, 0, 0, 0, 0, 0, 0],
49
           [0, 0, 0, 0, 0, 0, 0],
50
           [0, 0, 0, 0, 0, 0, 0]])
51
    >>> corner_moravec(square).astype(int)
52
    array([[0, 0, 0, 0, 0, 0, 0],
53
           [0, 0, 0, 0, 0, 0, 0],
54
           [0, 0, 1, 1, 1, 0, 0],
55
           [0, 0, 1, 2, 1, 0, 0],
56
           [0, 0, 1, 1, 1, 0, 0],
57
           [0, 0, 0, 0, 0, 0, 0],
58
           [0, 0, 0, 0, 0, 0, 0]])
59
    """
60

61
    cdef Py_ssize_t rows = cimage.shape[0]
62
    cdef Py_ssize_t cols = cimage.shape[1]
63

64
    if np_floats is cnp.float32_t:
65
        dtype = np.float32
66
    else:
67
        dtype = np.float64
68

69
    cdef np_floats[:, ::1] out = np.zeros((rows, cols), dtype=dtype)
70

71
    cdef np_floats msum, min_msum, t
72
    cdef Py_ssize_t r, c, br, bc, mr, mc
73

74
    with nogil:
75
        for r in range(2 * window_size, rows - 2 * window_size):
76
            for c in range(2 * window_size, cols - 2 * window_size):
77
                min_msum = DBL_MAX
78
                for br in range(r - window_size, r + window_size + 1):
79
                    for bc in range(c - window_size, c + window_size + 1):
80
                        if br != r and bc != c:
81
                            msum = 0
82
                            for mr in range(- window_size, window_size + 1):
83
                                for mc in range(- window_size, window_size + 1):
84
                                    t = cimage[r + mr, c + mc] - cimage[br + mr, bc + mc]
85
                                    msum += t * t
86
                            min_msum = min(msum, min_msum)
87

88
                out[r, c] = min_msum
89

90
    return np.asarray(out)
91

92

93
cdef inline np_floats _corner_fast_response(np_floats curr_pixel,
94
                                            np_floats* circle_intensities,
95
                                            signed char* bins, signed char
96
                                            state, char n) noexcept nogil:
97
    cdef char consecutive_count = 0
98
    cdef np_floats curr_response
99
    cdef Py_ssize_t l, m
100
    for l in range(15 + n):
101
        if bins[l % 16] == state:
102
            consecutive_count += 1
103
            if consecutive_count == n:
104
                curr_response = 0
105
                for m in range(16):
106
                    curr_response += fabs(circle_intensities[m] - curr_pixel)
107
                return curr_response
108
        else:
109
            consecutive_count = 0
110
    return 0
111

112

113
def _corner_fast(np_floats[:, ::1] image, signed char n, np_floats threshold):
114

115
    if np_floats is cnp.float32_t:
116
        dtype = np.float32
117
    else:
118
        dtype = np.float64
119

120
    cdef Py_ssize_t rows = image.shape[0]
121
    cdef Py_ssize_t cols = image.shape[1]
122

123
    cdef Py_ssize_t i, j, k
124

125
    cdef signed char speed_sum_b, speed_sum_d
126
    cdef np_floats curr_pixel, ring_pixel
127
    cdef np_floats lower_threshold, upper_threshold
128
    cdef np_floats[:, ::1] corner_response = np.zeros((rows, cols),
129
                                                      dtype=dtype)
130

131
    cdef signed char *rp = [0, 1, 2, 3, 3, 3, 2, 1, 0, -1, -2, -3, -3,
132
                            -3, -2, -1]
133
    cdef signed char *cp = [3, 3, 2, 1, 0, -1, -2, -3, -3, -3, -2, -1,
134
                            0, 1, 2, 3]
135
    cdef signed char bins[16]
136
    cdef np_floats circle_intensities[16]
137

138
    cdef cnp.float64_t curr_response
139

140
    with nogil:
141
        for i in range(3, rows - 3):
142
            for j in range(3, cols - 3):
143

144
                curr_pixel = image[i, j]
145
                lower_threshold = curr_pixel - threshold
146
                upper_threshold = curr_pixel + threshold
147

148
                # High speed test for n >= 12
149
                if n >= 12:
150
                    speed_sum_b = 0
151
                    speed_sum_d = 0
152
                    for k in range(0, 16, 4):
153
                        ring_pixel = image[i + rp[k], j + cp[k]]
154
                        if ring_pixel > upper_threshold:
155
                            speed_sum_b += 1
156
                        elif ring_pixel < lower_threshold:
157
                            speed_sum_d += 1
158
                    if speed_sum_d < 3 and speed_sum_b < 3:
159
                        continue
160

161
                for k in range(16):
162
                    circle_intensities[k] = image[i + rp[k], j + cp[k]]
163
                    if circle_intensities[k] > upper_threshold:
164
                        # Brighter pixel
165
                        bins[k] = b'b'
166
                    elif circle_intensities[k] < lower_threshold:
167
                        # Darker pixel
168
                        bins[k] = b'd'
169
                    else:
170
                        # Similar pixel
171
                        bins[k] = b's'
172

173
                # Test for bright pixels
174
                curr_response = _corner_fast_response[np_floats](curr_pixel,
175
                                                      circle_intensities, bins,
176
                                                      b'b', n)
177

178
                # Test for dark pixels
179
                if curr_response == 0:
180
                    curr_response = _corner_fast_response[np_floats](curr_pixel,
181
                                                          circle_intensities,
182
                                                          bins, b'd', n)
183

184
                corner_response[i, j] = curr_response
185

186
    return np.asarray(corner_response)
187

188

189
def _corner_orientations(np_floats[:, ::1] image, Py_ssize_t[:, :] corners,
190
                         mask):
191
    """Compute the orientation of corners.
192

193
    The orientation of corners is computed using the first order central moment
194
    i.e. the center of mass approach. The corner orientation is the angle of
195
    the vector from the corner coordinate to the intensity centroid in the
196
    local neighborhood around the corner calculated using first order central
197
    moment.
198

199
    Parameters
200
    ----------
201
    image : 2D array
202
        Input grayscale image.
203
    corners : (N, 2) array
204
        Corner coordinates as ``(row, col)``.
205
    mask : 2D array
206
        Mask defining the local neighborhood of the corner used for the
207
        calculation of the central moment.
208

209
    Returns
210
    -------
211
    orientations : (N, 1) array
212
        Orientations of corners in the range [-pi, pi].
213

214
    References
215
    ----------
216
    .. [1] Ethan Rublee, Vincent Rabaud, Kurt Konolige and Gary Bradski
217
          "ORB : An efficient alternative to SIFT and SURF"
218
          http://www.vision.cs.chubu.ac.jp/CV-R/pdf/Rublee_iccv2011.pdf
219
    .. [2] Paul L. Rosin, "Measuring Corner Properties"
220
          http://users.cs.cf.ac.uk/Paul.Rosin/corner2.pdf
221

222
    Examples
223
    --------
224
    >>> from skimage.morphology import octagon
225
    >>> from skimage.feature import (corner_fast, corner_peaks,
226
    ...                              corner_orientations)
227
    >>> square = np.zeros((12, 12))
228
    >>> square[3:9, 3:9] = 1
229
    >>> square.astype(int)
230
    array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
231
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
232
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
233
           [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0],
234
           [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0],
235
           [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0],
236
           [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0],
237
           [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0],
238
           [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0],
239
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
240
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
241
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
242
    >>> corners = corner_peaks(corner_fast(square, 9), min_distance=1)
243
    >>> corners
244
    array([[3, 3],
245
           [3, 8],
246
           [8, 3],
247
           [8, 8]])
248
    >>> orientations = corner_orientations(square, corners, octagon(3, 2))
249
    >>> np.rad2deg(orientations)
250
    array([  45.,  135.,  -45., -135.])
251

252
    """
253

254
    if np_floats is cnp.float32_t:
255
        dtype = np.float32
256
    else:
257
        dtype = np.float64
258

259
    if mask.shape[0] % 2 != 1 or mask.shape[1] % 2 != 1:
260
        raise ValueError("Size of mask must be uneven.")
261

262
    cdef unsigned char[:, ::1] cmask = np.ascontiguousarray(mask != 0,
263
                                                            dtype=np.uint8)
264

265
    cdef Py_ssize_t i, r, c, r0, c0
266
    cdef Py_ssize_t mrows = mask.shape[0]
267
    cdef Py_ssize_t mcols = mask.shape[1]
268
    cdef Py_ssize_t mrows2 = (mrows - 1) / 2
269
    cdef Py_ssize_t mcols2 = (mcols - 1) / 2
270
    cdef np_floats[:, :] cimage = np.pad(image, (mrows2, mcols2),
271
                                         mode='constant',
272
                                         constant_values=0)
273
    cdef np_floats[:] orientations = np.zeros(corners.shape[0], dtype=dtype)
274
    cdef np_floats curr_pixel, m01, m10, m01_tmp
275

276
    with nogil:
277
        for i in range(corners.shape[0]):
278
            r0 = corners[i, 0]
279
            c0 = corners[i, 1]
280

281
            m01 = 0
282
            m10 = 0
283

284
            for r in range(mrows):
285
                m01_tmp = 0
286
                for c in range(mcols):
287
                    if cmask[r, c]:
288
                        curr_pixel = cimage[r0 + r, c0 + c]
289
                        m10 += curr_pixel * (c - mcols2)
290
                        m01_tmp += curr_pixel
291
                m01 += m01_tmp * (r - mrows2)
292

293
            orientations[i] = atan2(m01, m10)
294

295
    return np.asarray(orientations)
296