scikit-image

1
import numpy as np
2
from scipy.ndimage import maximum_filter, minimum_filter, convolve
3

4
from ..transform import integral_image
5
from .corner import structure_tensor
6
from ..morphology import octagon, star
7
from .censure_cy import _censure_dob_loop
8
from ..feature.util import (
9
    FeatureDetector,
10
    _prepare_grayscale_input_2D,
11
    _mask_border_keypoints,
12
)
13
from .._shared.utils import check_nD
14

15
# The paper(Reference [1]) mentions the sizes of the Octagon shaped filter
16
# kernel for the first seven scales only. The sizes of the later scales
17
# have been extrapolated based on the following statement in the paper.
18
# "These octagons scale linearly and were experimentally chosen to correspond
19
# to the seven DOBs described in the previous section."
20
OCTAGON_OUTER_SHAPE = [
21
    (5, 2),
22
    (5, 3),
23
    (7, 3),
24
    (9, 4),
25
    (9, 7),
26
    (13, 7),
27
    (15, 10),
28
    (15, 11),
29
    (15, 12),
30
    (17, 13),
31
    (17, 14),
32
]
33
OCTAGON_INNER_SHAPE = [
34
    (3, 0),
35
    (3, 1),
36
    (3, 2),
37
    (5, 2),
38
    (5, 3),
39
    (5, 4),
40
    (5, 5),
41
    (7, 5),
42
    (7, 6),
43
    (9, 6),
44
    (9, 7),
45
]
46

47
# The sizes for the STAR shaped filter kernel for different scales have been
48
# taken from the OpenCV implementation.
49
STAR_SHAPE = [1, 2, 3, 4, 6, 8, 11, 12, 16, 22, 23, 32, 45, 46, 64, 90, 128]
50
STAR_FILTER_SHAPE = [
51
    (1, 0),
52
    (3, 1),
53
    (4, 2),
54
    (5, 3),
55
    (7, 4),
56
    (8, 5),
57
    (9, 6),
58
    (11, 8),
59
    (13, 10),
60
    (14, 11),
61
    (15, 12),
62
    (16, 14),
63
]
64

65

66
def _filter_image(image, min_scale, max_scale, mode):
67
    response = np.zeros(
68
        (image.shape[0], image.shape[1], max_scale - min_scale + 1), dtype=np.float64
69
    )
70

71
    if mode == 'dob':
72
        # make response[:, :, i] contiguous memory block
73
        item_size = response.itemsize
74
        response.strides = (
75
            item_size * response.shape[1],
76
            item_size,
77
            item_size * response.shape[0] * response.shape[1],
78
        )
79

80
        integral_img = integral_image(image)
81

82
        for i in range(max_scale - min_scale + 1):
83
            n = min_scale + i
84

85
            # Constant multipliers for the outer region and the inner region
86
            # of the bi-level filters with the constraint of keeping the
87
            # DC bias 0.
88
            inner_weight = 1.0 / (2 * n + 1) ** 2
89
            outer_weight = 1.0 / (12 * n**2 + 4 * n)
90

91
            _censure_dob_loop(
92
                n, integral_img, response[:, :, i], inner_weight, outer_weight
93
            )
94

95
    # NOTE : For the Octagon shaped filter, we implemented and evaluated the
96
    # slanted integral image based image filtering but the performance was
97
    # more or less equal to image filtering using
98
    # scipy.ndimage.filters.convolve(). Hence we have decided to use the
99
    # later for a much cleaner implementation.
100
    elif mode == 'octagon':
101
        # TODO : Decide the shapes of Octagon filters for scales > 7
102

103
        for i in range(max_scale - min_scale + 1):
104
            mo, no = OCTAGON_OUTER_SHAPE[min_scale + i - 1]
105
            mi, ni = OCTAGON_INNER_SHAPE[min_scale + i - 1]
106
            response[:, :, i] = convolve(image, _octagon_kernel(mo, no, mi, ni))
107

108
    elif mode == 'star':
109
        for i in range(max_scale - min_scale + 1):
110
            m = STAR_SHAPE[STAR_FILTER_SHAPE[min_scale + i - 1][0]]
111
            n = STAR_SHAPE[STAR_FILTER_SHAPE[min_scale + i - 1][1]]
112
            response[:, :, i] = convolve(image, _star_kernel(m, n))
113

114
    return response
115

116

117
def _octagon_kernel(mo, no, mi, ni):
118
    outer = (mo + 2 * no) ** 2 - 2 * no * (no + 1)
119
    inner = (mi + 2 * ni) ** 2 - 2 * ni * (ni + 1)
120
    outer_weight = 1.0 / (outer - inner)
121
    inner_weight = 1.0 / inner
122
    c = ((mo + 2 * no) - (mi + 2 * ni)) // 2
123
    outer_oct = octagon(mo, no)
124
    inner_oct = np.zeros((mo + 2 * no, mo + 2 * no))
125
    inner_oct[c:-c, c:-c] = octagon(mi, ni)
126
    bfilter = outer_weight * outer_oct - (outer_weight + inner_weight) * inner_oct
127
    return bfilter
128

129

130
def _star_kernel(m, n):
131
    c = m + m // 2 - n - n // 2
132
    outer_star = star(m)
133
    inner_star = np.zeros_like(outer_star)
134
    inner_star[c:-c, c:-c] = star(n)
135
    outer_weight = 1.0 / (np.sum(outer_star - inner_star))
136
    inner_weight = 1.0 / np.sum(inner_star)
137
    bfilter = outer_weight * outer_star - (outer_weight + inner_weight) * inner_star
138
    return bfilter
139

140

141
def _suppress_lines(feature_mask, image, sigma, line_threshold):
142
    Arr, Arc, Acc = structure_tensor(image, sigma, order='rc')
143
    feature_mask[(Arr + Acc) ** 2 > line_threshold * (Arr * Acc - Arc**2)] = False
144

145

146
class CENSURE(FeatureDetector):
147
    """CENSURE keypoint detector.
148

149
    min_scale : int, optional
150
        Minimum scale to extract keypoints from.
151
    max_scale : int, optional
152
        Maximum scale to extract keypoints from. The keypoints will be
153
        extracted from all the scales except the first and the last i.e.
154
        from the scales in the range [min_scale + 1, max_scale - 1]. The filter
155
        sizes for different scales is such that the two adjacent scales
156
        comprise of an octave.
157
    mode : {'DoB', 'Octagon', 'STAR'}, optional
158
        Type of bi-level filter used to get the scales of the input image.
159
        Possible values are 'DoB', 'Octagon' and 'STAR'. The three modes
160
        represent the shape of the bi-level filters i.e. box(square), octagon
161
        and star respectively. For instance, a bi-level octagon filter consists
162
        of a smaller inner octagon and a larger outer octagon with the filter
163
        weights being uniformly negative in both the inner octagon while
164
        uniformly positive in the difference region. Use STAR and Octagon for
165
        better features and DoB for better performance.
166
    non_max_threshold : float, optional
167
        Threshold value used to suppress maximas and minimas with a weak
168
        magnitude response obtained after Non-Maximal Suppression.
169
    line_threshold : float, optional
170
        Threshold for rejecting interest points which have ratio of principal
171
        curvatures greater than this value.
172

173
    Attributes
174
    ----------
175
    keypoints : (N, 2) array
176
        Keypoint coordinates as ``(row, col)``.
177
    scales : (N,) array
178
        Corresponding scales.
179

180
    References
181
    ----------
182
    .. [1] Motilal Agrawal, Kurt Konolige and Morten Rufus Blas
183
           "CENSURE: Center Surround Extremas for Realtime Feature
184
           Detection and Matching",
185
           https://link.springer.com/chapter/10.1007/978-3-540-88693-8_8
186
           :DOI:`10.1007/978-3-540-88693-8_8`
187

188
    .. [2] Adam Schmidt, Marek Kraft, Michal Fularz and Zuzanna Domagala
189
           "Comparative Assessment of Point Feature Detectors and
190
           Descriptors in the Context of Robot Navigation"
191
           http://yadda.icm.edu.pl/yadda/element/bwmeta1.element.baztech-268aaf28-0faf-4872-a4df-7e2e61cb364c/c/Schmidt_comparative.pdf
192
           :DOI:`10.1.1.465.1117`
193

194
    Examples
195
    --------
196
    >>> from skimage.data import astronaut
197
    >>> from skimage.color import rgb2gray
198
    >>> from skimage.feature import CENSURE
199
    >>> img = rgb2gray(astronaut()[100:300, 100:300])
200
    >>> censure = CENSURE()
201
    >>> censure.detect(img)
202
    >>> censure.keypoints
203
    array([[  4, 148],
204
           [ 12,  73],
205
           [ 21, 176],
206
           [ 91,  22],
207
           [ 93,  56],
208
           [ 94,  22],
209
           [ 95,  54],
210
           [100,  51],
211
           [103,  51],
212
           [106,  67],
213
           [108,  15],
214
           [117,  20],
215
           [122,  60],
216
           [125,  37],
217
           [129,  37],
218
           [133,  76],
219
           [145,  44],
220
           [146,  94],
221
           [150, 114],
222
           [153,  33],
223
           [154, 156],
224
           [155, 151],
225
           [184,  63]])
226
    >>> censure.scales
227
    array([2, 6, 6, 2, 4, 3, 2, 3, 2, 6, 3, 2, 2, 3, 2, 2, 2, 3, 2, 2, 4, 2,
228
           2])
229

230
    """
231

232
    def __init__(
233
        self,
234
        min_scale=1,
235
        max_scale=7,
236
        mode='DoB',
237
        non_max_threshold=0.15,
238
        line_threshold=10,
239
    ):
240
        mode = mode.lower()
241
        if mode not in ('dob', 'octagon', 'star'):
242
            raise ValueError("`mode` must be one of 'DoB', 'Octagon', 'STAR'.")
243

244
        if min_scale < 1 or max_scale < 1 or max_scale - min_scale < 2:
245
            raise ValueError(
246
                'The scales must be >= 1 and the number of ' 'scales should be >= 3.'
247
            )
248

249
        self.min_scale = min_scale
250
        self.max_scale = max_scale
251
        self.mode = mode
252
        self.non_max_threshold = non_max_threshold
253
        self.line_threshold = line_threshold
254

255
        self.keypoints = None
256
        self.scales = None
257

258
    def detect(self, image):
259
        """Detect CENSURE keypoints along with the corresponding scale.
260

261
        Parameters
262
        ----------
263
        image : 2D ndarray
264
            Input image.
265

266
        """
267

268
        # (1) First we generate the required scales on the input grayscale
269
        # image using a bi-level filter and stack them up in `filter_response`.
270

271
        # (2) We then perform Non-Maximal suppression in 3 x 3 x 3 window on
272
        # the filter_response to suppress points that are neither minima or
273
        # maxima in 3 x 3 x 3 neighborhood. We obtain a boolean ndarray
274
        # `feature_mask` containing all the minimas and maximas in
275
        # `filter_response` as True.
276
        # (3) Then we suppress all the points in the `feature_mask` for which
277
        # the corresponding point in the image at a particular scale has the
278
        # ratio of principal curvatures greater than `line_threshold`.
279
        # (4) Finally, we remove the border keypoints and return the keypoints
280
        # along with its corresponding scale.
281

282
        check_nD(image, 2)
283

284
        num_scales = self.max_scale - self.min_scale
285

286
        image = np.ascontiguousarray(_prepare_grayscale_input_2D(image))
287

288
        # Generating all the scales
289
        filter_response = _filter_image(
290
            image, self.min_scale, self.max_scale, self.mode
291
        )
292

293
        # Suppressing points that are neither minima or maxima in their
294
        # 3 x 3 x 3 neighborhood to zero
295
        minimas = minimum_filter(filter_response, (3, 3, 3)) == filter_response
296
        maximas = maximum_filter(filter_response, (3, 3, 3)) == filter_response
297

298
        feature_mask = minimas | maximas
299
        feature_mask[filter_response < self.non_max_threshold] = False
300

301
        for i in range(1, num_scales):
302
            # sigma = (window_size - 1) / 6.0, so the window covers > 99% of
303
            #                                  the kernel's distribution
304
            # window_size = 7 + 2 * (min_scale - 1 + i)
305
            # Hence sigma = 1 + (min_scale - 1 + i)/ 3.0
306
            _suppress_lines(
307
                feature_mask[:, :, i],
308
                image,
309
                (1 + (self.min_scale + i - 1) / 3.0),
310
                self.line_threshold,
311
            )
312

313
        rows, cols, scales = np.nonzero(feature_mask[..., 1:num_scales])
314
        keypoints = np.column_stack([rows, cols])
315
        scales = scales + self.min_scale + 1
316

317
        if self.mode == 'dob':
318
            self.keypoints = keypoints
319
            self.scales = scales
320
            return
321

322
        cumulative_mask = np.zeros(keypoints.shape[0], dtype=bool)
323

324
        if self.mode == 'octagon':
325
            for i in range(self.min_scale + 1, self.max_scale):
326
                c = (OCTAGON_OUTER_SHAPE[i - 1][0] - 1) // 2 + OCTAGON_OUTER_SHAPE[
327
                    i - 1
328
                ][1]
329
                cumulative_mask |= _mask_border_keypoints(image.shape, keypoints, c) & (
330
                    scales == i
331
                )
332
        elif self.mode == 'star':
333
            for i in range(self.min_scale + 1, self.max_scale):
334
                c = (
335
                    STAR_SHAPE[STAR_FILTER_SHAPE[i - 1][0]]
336
                    + STAR_SHAPE[STAR_FILTER_SHAPE[i - 1][0]] // 2
337
                )
338
                cumulative_mask |= _mask_border_keypoints(image.shape, keypoints, c) & (
339
                    scales == i
340
                )
341

342
        self.keypoints = keypoints[cumulative_mask]
343
        self.scales = scales[cumulative_mask]
344