censure.py 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346
  1. import numpy as np
  2. from scipy.ndimage import maximum_filter, minimum_filter, convolve
  3. from ..transform import integral_image
  4. from .corner import structure_tensor
  5. from ..morphology import octagon, star
  6. from .censure_cy import _censure_dob_loop
  7. from ..feature.util import (
  8. FeatureDetector,
  9. _prepare_grayscale_input_2D,
  10. _mask_border_keypoints,
  11. )
  12. from .._shared.utils import check_nD
  13. # The paper(Reference [1]) mentions the sizes of the Octagon shaped filter
  14. # kernel for the first seven scales only. The sizes of the later scales
  15. # have been extrapolated based on the following statement in the paper.
  16. # "These octagons scale linearly and were experimentally chosen to correspond
  17. # to the seven DOBs described in the previous section."
  18. OCTAGON_OUTER_SHAPE = [
  19. (5, 2),
  20. (5, 3),
  21. (7, 3),
  22. (9, 4),
  23. (9, 7),
  24. (13, 7),
  25. (15, 10),
  26. (15, 11),
  27. (15, 12),
  28. (17, 13),
  29. (17, 14),
  30. ]
  31. OCTAGON_INNER_SHAPE = [
  32. (3, 0),
  33. (3, 1),
  34. (3, 2),
  35. (5, 2),
  36. (5, 3),
  37. (5, 4),
  38. (5, 5),
  39. (7, 5),
  40. (7, 6),
  41. (9, 6),
  42. (9, 7),
  43. ]
  44. # The sizes for the STAR shaped filter kernel for different scales have been
  45. # taken from the OpenCV implementation.
  46. STAR_SHAPE = [1, 2, 3, 4, 6, 8, 11, 12, 16, 22, 23, 32, 45, 46, 64, 90, 128]
  47. STAR_FILTER_SHAPE = [
  48. (1, 0),
  49. (3, 1),
  50. (4, 2),
  51. (5, 3),
  52. (7, 4),
  53. (8, 5),
  54. (9, 6),
  55. (11, 8),
  56. (13, 10),
  57. (14, 11),
  58. (15, 12),
  59. (16, 14),
  60. ]
  61. def _filter_image(image, min_scale, max_scale, mode):
  62. response = np.zeros(
  63. (image.shape[0], image.shape[1], max_scale - min_scale + 1), dtype=np.float64
  64. )
  65. if mode == 'dob':
  66. # make response[:, :, i] contiguous memory block
  67. item_size = response.itemsize
  68. response = np.lib.stride_tricks.as_strided(
  69. response,
  70. strides=(
  71. item_size * response.shape[1],
  72. item_size,
  73. item_size * response.shape[0] * response.shape[1],
  74. ),
  75. )
  76. integral_img = integral_image(image)
  77. for i in range(max_scale - min_scale + 1):
  78. n = min_scale + i
  79. # Constant multipliers for the outer region and the inner region
  80. # of the bi-level filters with the constraint of keeping the
  81. # DC bias 0.
  82. inner_weight = 1.0 / (2 * n + 1) ** 2
  83. outer_weight = 1.0 / (12 * n**2 + 4 * n)
  84. _censure_dob_loop(
  85. n, integral_img, response[:, :, i], inner_weight, outer_weight
  86. )
  87. # NOTE : For the Octagon shaped filter, we implemented and evaluated the
  88. # slanted integral image based image filtering but the performance was
  89. # more or less equal to image filtering using
  90. # scipy.ndimage.filters.convolve(). Hence we have decided to use the
  91. # later for a much cleaner implementation.
  92. elif mode == 'octagon':
  93. # TODO : Decide the shapes of Octagon filters for scales > 7
  94. for i in range(max_scale - min_scale + 1):
  95. mo, no = OCTAGON_OUTER_SHAPE[min_scale + i - 1]
  96. mi, ni = OCTAGON_INNER_SHAPE[min_scale + i - 1]
  97. response[:, :, i] = convolve(image, _octagon_kernel(mo, no, mi, ni))
  98. elif mode == 'star':
  99. for i in range(max_scale - min_scale + 1):
  100. m = STAR_SHAPE[STAR_FILTER_SHAPE[min_scale + i - 1][0]]
  101. n = STAR_SHAPE[STAR_FILTER_SHAPE[min_scale + i - 1][1]]
  102. response[:, :, i] = convolve(image, _star_kernel(m, n))
  103. return response
  104. def _octagon_kernel(mo, no, mi, ni):
  105. outer = (mo + 2 * no) ** 2 - 2 * no * (no + 1)
  106. inner = (mi + 2 * ni) ** 2 - 2 * ni * (ni + 1)
  107. outer_weight = 1.0 / (outer - inner)
  108. inner_weight = 1.0 / inner
  109. c = ((mo + 2 * no) - (mi + 2 * ni)) // 2
  110. outer_oct = octagon(mo, no)
  111. inner_oct = np.zeros((mo + 2 * no, mo + 2 * no))
  112. inner_oct[c:-c, c:-c] = octagon(mi, ni)
  113. bfilter = outer_weight * outer_oct - (outer_weight + inner_weight) * inner_oct
  114. return bfilter
  115. def _star_kernel(m, n):
  116. c = m + m // 2 - n - n // 2
  117. outer_star = star(m)
  118. inner_star = np.zeros_like(outer_star)
  119. inner_star[c:-c, c:-c] = star(n)
  120. outer_weight = 1.0 / (np.sum(outer_star - inner_star))
  121. inner_weight = 1.0 / np.sum(inner_star)
  122. bfilter = outer_weight * outer_star - (outer_weight + inner_weight) * inner_star
  123. return bfilter
  124. def _suppress_lines(feature_mask, image, sigma, line_threshold):
  125. Arr, Arc, Acc = structure_tensor(image, sigma, order='rc')
  126. feature_mask[(Arr + Acc) ** 2 > line_threshold * (Arr * Acc - Arc**2)] = False
  127. class CENSURE(FeatureDetector):
  128. """CENSURE keypoint detector.
  129. min_scale : int, optional
  130. Minimum scale to extract keypoints from.
  131. max_scale : int, optional
  132. Maximum scale to extract keypoints from. The keypoints will be
  133. extracted from all the scales except the first and the last i.e.
  134. from the scales in the range [min_scale + 1, max_scale - 1]. The filter
  135. sizes for different scales is such that the two adjacent scales
  136. comprise of an octave.
  137. mode : {'DoB', 'Octagon', 'STAR'}, optional
  138. Type of bi-level filter used to get the scales of the input image.
  139. Possible values are 'DoB', 'Octagon' and 'STAR'. The three modes
  140. represent the shape of the bi-level filters i.e. box(square), octagon
  141. and star respectively. For instance, a bi-level octagon filter consists
  142. of a smaller inner octagon and a larger outer octagon with the filter
  143. weights being uniformly negative in both the inner octagon while
  144. uniformly positive in the difference region. Use STAR and Octagon for
  145. better features and DoB for better performance.
  146. non_max_threshold : float, optional
  147. Threshold value used to suppress maximas and minimas with a weak
  148. magnitude response obtained after Non-Maximal Suppression.
  149. line_threshold : float, optional
  150. Threshold for rejecting interest points which have ratio of principal
  151. curvatures greater than this value.
  152. Attributes
  153. ----------
  154. keypoints : (N, 2) array
  155. Keypoint coordinates as ``(row, col)``.
  156. scales : (N,) array
  157. Corresponding scales.
  158. References
  159. ----------
  160. .. [1] Motilal Agrawal, Kurt Konolige and Morten Rufus Blas
  161. "CENSURE: Center Surround Extremas for Realtime Feature
  162. Detection and Matching",
  163. https://link.springer.com/chapter/10.1007/978-3-540-88693-8_8
  164. :DOI:`10.1007/978-3-540-88693-8_8`
  165. .. [2] Adam Schmidt, Marek Kraft, Michal Fularz and Zuzanna Domagala
  166. "Comparative Assessment of Point Feature Detectors and
  167. Descriptors in the Context of Robot Navigation"
  168. http://yadda.icm.edu.pl/yadda/element/bwmeta1.element.baztech-268aaf28-0faf-4872-a4df-7e2e61cb364c/c/Schmidt_comparative.pdf
  169. :DOI:`10.1.1.465.1117`
  170. Examples
  171. --------
  172. >>> from skimage.data import astronaut
  173. >>> from skimage.color import rgb2gray
  174. >>> from skimage.feature import CENSURE
  175. >>> img = rgb2gray(astronaut()[100:300, 100:300])
  176. >>> censure = CENSURE()
  177. >>> censure.detect(img)
  178. >>> censure.keypoints
  179. array([[ 4, 148],
  180. [ 12, 73],
  181. [ 21, 176],
  182. [ 91, 22],
  183. [ 93, 56],
  184. [ 94, 22],
  185. [ 95, 54],
  186. [100, 51],
  187. [103, 51],
  188. [106, 67],
  189. [108, 15],
  190. [117, 20],
  191. [122, 60],
  192. [125, 37],
  193. [129, 37],
  194. [133, 76],
  195. [145, 44],
  196. [146, 94],
  197. [150, 114],
  198. [153, 33],
  199. [154, 156],
  200. [155, 151],
  201. [184, 63]])
  202. >>> censure.scales
  203. array([2, 6, 6, 2, 4, 3, 2, 3, 2, 6, 3, 2, 2, 3, 2, 2, 2, 3, 2, 2, 4, 2,
  204. 2])
  205. """
  206. def __init__(
  207. self,
  208. min_scale=1,
  209. max_scale=7,
  210. mode='DoB',
  211. non_max_threshold=0.15,
  212. line_threshold=10,
  213. ):
  214. mode = mode.lower()
  215. if mode not in ('dob', 'octagon', 'star'):
  216. raise ValueError("`mode` must be one of 'DoB', 'Octagon', 'STAR'.")
  217. if min_scale < 1 or max_scale < 1 or max_scale - min_scale < 2:
  218. raise ValueError(
  219. 'The scales must be >= 1 and the number of ' 'scales should be >= 3.'
  220. )
  221. self.min_scale = min_scale
  222. self.max_scale = max_scale
  223. self.mode = mode
  224. self.non_max_threshold = non_max_threshold
  225. self.line_threshold = line_threshold
  226. self.keypoints = None
  227. self.scales = None
  228. def detect(self, image):
  229. """Detect CENSURE keypoints along with the corresponding scale.
  230. Parameters
  231. ----------
  232. image : 2D ndarray
  233. Input image.
  234. """
  235. # (1) First we generate the required scales on the input grayscale
  236. # image using a bi-level filter and stack them up in `filter_response`.
  237. # (2) We then perform Non-Maximal suppression in 3 x 3 x 3 window on
  238. # the filter_response to suppress points that are neither minima or
  239. # maxima in 3 x 3 x 3 neighborhood. We obtain a boolean ndarray
  240. # `feature_mask` containing all the minimas and maximas in
  241. # `filter_response` as True.
  242. # (3) Then we suppress all the points in the `feature_mask` for which
  243. # the corresponding point in the image at a particular scale has the
  244. # ratio of principal curvatures greater than `line_threshold`.
  245. # (4) Finally, we remove the border keypoints and return the keypoints
  246. # along with its corresponding scale.
  247. check_nD(image, 2)
  248. num_scales = self.max_scale - self.min_scale
  249. image = np.ascontiguousarray(_prepare_grayscale_input_2D(image))
  250. # Generating all the scales
  251. filter_response = _filter_image(
  252. image, self.min_scale, self.max_scale, self.mode
  253. )
  254. # Suppressing points that are neither minima or maxima in their
  255. # 3 x 3 x 3 neighborhood to zero
  256. minimas = minimum_filter(filter_response, (3, 3, 3)) == filter_response
  257. maximas = maximum_filter(filter_response, (3, 3, 3)) == filter_response
  258. feature_mask = minimas | maximas
  259. feature_mask[filter_response < self.non_max_threshold] = False
  260. for i in range(1, num_scales):
  261. # sigma = (window_size - 1) / 6.0, so the window covers > 99% of
  262. # the kernel's distribution
  263. # window_size = 7 + 2 * (min_scale - 1 + i)
  264. # Hence sigma = 1 + (min_scale - 1 + i)/ 3.0
  265. _suppress_lines(
  266. feature_mask[:, :, i],
  267. image,
  268. (1 + (self.min_scale + i - 1) / 3.0),
  269. self.line_threshold,
  270. )
  271. rows, cols, scales = np.nonzero(feature_mask[..., 1:num_scales])
  272. keypoints = np.column_stack([rows, cols])
  273. scales = scales + self.min_scale + 1
  274. if self.mode == 'dob':
  275. self.keypoints = keypoints
  276. self.scales = scales
  277. return
  278. cumulative_mask = np.zeros(keypoints.shape[0], dtype=bool)
  279. if self.mode == 'octagon':
  280. for i in range(self.min_scale + 1, self.max_scale):
  281. c = (OCTAGON_OUTER_SHAPE[i - 1][0] - 1) // 2 + OCTAGON_OUTER_SHAPE[
  282. i - 1
  283. ][1]
  284. cumulative_mask |= _mask_border_keypoints(image.shape, keypoints, c) & (
  285. scales == i
  286. )
  287. elif self.mode == 'star':
  288. for i in range(self.min_scale + 1, self.max_scale):
  289. c = (
  290. STAR_SHAPE[STAR_FILTER_SHAPE[i - 1][0]]
  291. + STAR_SHAPE[STAR_FILTER_SHAPE[i - 1][0]] // 2
  292. )
  293. cumulative_mask |= _mask_border_keypoints(image.shape, keypoints, c) & (
  294. scales == i
  295. )
  296. self.keypoints = keypoints[cumulative_mask]
  297. self.scales = scales[cumulative_mask]