_daisy.py 9.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249
  1. import math
  2. import numpy as np
  3. from numpy import arctan2, exp, pi, sqrt
  4. from .. import draw
  5. from ..util.dtype import img_as_float
  6. from .._shared.filters import gaussian
  7. from .._shared.utils import check_nD
  8. from ..color import gray2rgb
  9. def daisy(
  10. image,
  11. step=4,
  12. radius=15,
  13. rings=3,
  14. histograms=8,
  15. orientations=8,
  16. normalization='l1',
  17. sigmas=None,
  18. ring_radii=None,
  19. visualize=False,
  20. ):
  21. '''Extract DAISY feature descriptors densely for the given image.
  22. DAISY is a feature descriptor similar to SIFT formulated in a way that
  23. allows for fast dense extraction. Typically, this is practical for
  24. bag-of-features image representations.
  25. The implementation follows Tola et al. [1]_ but deviate on the following
  26. points:
  27. * Histogram bin contribution are smoothed with a circular Gaussian
  28. window over the tonal range (the angular range).
  29. * The sigma values of the spatial Gaussian smoothing in this code do not
  30. match the sigma values in the original code by Tola et al. [2]_. In
  31. their code, spatial smoothing is applied to both the input image and
  32. the center histogram. However, this smoothing is not documented in [1]_
  33. and, therefore, it is omitted.
  34. Parameters
  35. ----------
  36. image : (M, N) array
  37. Input image (grayscale).
  38. step : int, optional
  39. Distance between descriptor sampling points.
  40. radius : int, optional
  41. Radius (in pixels) of the outermost ring.
  42. rings : int, optional
  43. Number of rings.
  44. histograms : int, optional
  45. Number of histograms sampled per ring.
  46. orientations : int, optional
  47. Number of orientations (bins) per histogram.
  48. normalization : [ 'l1' | 'l2' | 'daisy' | 'off' ], optional
  49. How to normalize the descriptors
  50. * 'l1': L1-normalization of each descriptor.
  51. * 'l2': L2-normalization of each descriptor.
  52. * 'daisy': L2-normalization of individual histograms.
  53. * 'off': Disable normalization.
  54. sigmas : 1D array of float, optional
  55. Standard deviation of spatial Gaussian smoothing for the center
  56. histogram and for each ring of histograms. The array of sigmas should
  57. be sorted from the center and out. I.e. the first sigma value defines
  58. the spatial smoothing of the center histogram and the last sigma value
  59. defines the spatial smoothing of the outermost ring. Specifying sigmas
  60. overrides the following parameter.
  61. ``rings = len(sigmas) - 1``
  62. ring_radii : 1D array of int, optional
  63. Radius (in pixels) for each ring. Specifying ring_radii overrides the
  64. following two parameters.
  65. ``rings = len(ring_radii)``
  66. ``radius = ring_radii[-1]``
  67. If both sigmas and ring_radii are given, they must satisfy the
  68. following predicate since no radius is needed for the center
  69. histogram.
  70. ``len(ring_radii) == len(sigmas) + 1``
  71. visualize : bool, optional
  72. Generate a visualization of the DAISY descriptors
  73. Returns
  74. -------
  75. descs : array
  76. Grid of DAISY descriptors for the given image as an array
  77. dimensionality (P, Q, R) where
  78. ``P = ceil((M - radius*2) / step)``
  79. ``Q = ceil((N - radius*2) / step)``
  80. ``R = (rings * histograms + 1) * orientations``
  81. descs_img : (M, N, 3) array (only if visualize==True)
  82. Visualization of the DAISY descriptors.
  83. References
  84. ----------
  85. .. [1] Tola et al. "Daisy: An efficient dense descriptor applied to wide-
  86. baseline stereo." Pattern Analysis and Machine Intelligence, IEEE
  87. Transactions on 32.5 (2010): 815-830.
  88. .. [2] http://cvlab.epfl.ch/software/daisy
  89. '''
  90. check_nD(image, 2, 'img')
  91. image = img_as_float(image)
  92. float_dtype = image.dtype
  93. # Validate parameters.
  94. if (
  95. sigmas is not None
  96. and ring_radii is not None
  97. and len(sigmas) - 1 != len(ring_radii)
  98. ):
  99. raise ValueError('`len(sigmas)-1 != len(ring_radii)`')
  100. if ring_radii is not None:
  101. rings = len(ring_radii)
  102. radius = ring_radii[-1]
  103. if sigmas is not None:
  104. rings = len(sigmas) - 1
  105. if sigmas is None:
  106. sigmas = [radius * (i + 1) / float(2 * rings) for i in range(rings)]
  107. if ring_radii is None:
  108. ring_radii = [radius * (i + 1) / float(rings) for i in range(rings)]
  109. if normalization not in ['l1', 'l2', 'daisy', 'off']:
  110. raise ValueError('Invalid normalization method.')
  111. # Compute image derivatives.
  112. dx = np.zeros(image.shape, dtype=float_dtype)
  113. dy = np.zeros(image.shape, dtype=float_dtype)
  114. dx[:, :-1] = np.diff(image, n=1, axis=1)
  115. dy[:-1, :] = np.diff(image, n=1, axis=0)
  116. # Compute gradient orientation and magnitude and their contribution
  117. # to the histograms.
  118. grad_mag = sqrt(dx**2 + dy**2)
  119. grad_ori = arctan2(dy, dx)
  120. orientation_kappa = orientations / pi
  121. orientation_angles = [2 * o * pi / orientations - pi for o in range(orientations)]
  122. hist = np.empty((orientations,) + image.shape, dtype=float_dtype)
  123. for i, o in enumerate(orientation_angles):
  124. # Weigh bin contribution by the circular normal distribution
  125. hist[i, :, :] = exp(orientation_kappa * np.cos(grad_ori - o))
  126. # Weigh bin contribution by the gradient magnitude
  127. hist[i, :, :] = np.multiply(hist[i, :, :], grad_mag)
  128. # Smooth orientation histograms for the center and all rings.
  129. sigmas = [sigmas[0]] + sigmas
  130. hist_smooth = np.empty((rings + 1,) + hist.shape, dtype=float_dtype)
  131. for i in range(rings + 1):
  132. for j in range(orientations):
  133. hist_smooth[i, j, :, :] = gaussian(
  134. hist[j, :, :], sigma=sigmas[i], mode='reflect'
  135. )
  136. # Assemble descriptor grid.
  137. theta = [2 * pi * j / histograms for j in range(histograms)]
  138. desc_dims = (rings * histograms + 1) * orientations
  139. descs = np.empty(
  140. (desc_dims, image.shape[0] - 2 * radius, image.shape[1] - 2 * radius),
  141. dtype=float_dtype,
  142. )
  143. descs[:orientations, :, :] = hist_smooth[0, :, radius:-radius, radius:-radius]
  144. idx = orientations
  145. for i in range(rings):
  146. for j in range(histograms):
  147. y_min = radius + int(round(ring_radii[i] * math.sin(theta[j])))
  148. y_max = descs.shape[1] + y_min
  149. x_min = radius + int(round(ring_radii[i] * math.cos(theta[j])))
  150. x_max = descs.shape[2] + x_min
  151. descs[idx : idx + orientations, :, :] = hist_smooth[
  152. i + 1, :, y_min:y_max, x_min:x_max
  153. ]
  154. idx += orientations
  155. descs = descs[:, ::step, ::step]
  156. descs = descs.swapaxes(0, 1).swapaxes(1, 2)
  157. # Normalize descriptors.
  158. if normalization != 'off':
  159. descs += 1e-10
  160. if normalization == 'l1':
  161. descs /= np.sum(descs, axis=2)[:, :, np.newaxis]
  162. elif normalization == 'l2':
  163. descs /= sqrt(np.sum(descs**2, axis=2))[:, :, np.newaxis]
  164. elif normalization == 'daisy':
  165. for i in range(0, desc_dims, orientations):
  166. norms = sqrt(np.sum(descs[:, :, i : i + orientations] ** 2, axis=2))
  167. descs[:, :, i : i + orientations] /= norms[:, :, np.newaxis]
  168. if visualize:
  169. descs_img = gray2rgb(image)
  170. for i in range(descs.shape[0]):
  171. for j in range(descs.shape[1]):
  172. # Draw center histogram sigma
  173. color = [1, 0, 0]
  174. desc_y = i * step + radius
  175. desc_x = j * step + radius
  176. rows, cols, val = draw.circle_perimeter_aa(
  177. desc_y, desc_x, int(sigmas[0])
  178. )
  179. draw.set_color(descs_img, (rows, cols), color, alpha=val)
  180. max_bin = np.max(descs[i, j, :])
  181. for o_num, o in enumerate(orientation_angles):
  182. # Draw center histogram bins
  183. bin_size = descs[i, j, o_num] / max_bin
  184. dy = sigmas[0] * bin_size * math.sin(o)
  185. dx = sigmas[0] * bin_size * math.cos(o)
  186. rows, cols, val = draw.line_aa(
  187. desc_y, desc_x, int(desc_y + dy), int(desc_x + dx)
  188. )
  189. draw.set_color(descs_img, (rows, cols), color, alpha=val)
  190. for r_num, r in enumerate(ring_radii):
  191. color_offset = float(1 + r_num) / rings
  192. color = (1 - color_offset, 1, color_offset)
  193. for t_num, t in enumerate(theta):
  194. # Draw ring histogram sigmas
  195. hist_y = desc_y + int(round(r * math.sin(t)))
  196. hist_x = desc_x + int(round(r * math.cos(t)))
  197. rows, cols, val = draw.circle_perimeter_aa(
  198. hist_y, hist_x, int(sigmas[r_num + 1])
  199. )
  200. draw.set_color(descs_img, (rows, cols), color, alpha=val)
  201. for o_num, o in enumerate(orientation_angles):
  202. # Draw histogram bins
  203. bin_size = descs[
  204. i,
  205. j,
  206. orientations
  207. + r_num * histograms * orientations
  208. + t_num * orientations
  209. + o_num,
  210. ]
  211. bin_size /= max_bin
  212. dy = sigmas[r_num + 1] * bin_size * math.sin(o)
  213. dx = sigmas[r_num + 1] * bin_size * math.cos(o)
  214. rows, cols, val = draw.line_aa(
  215. hist_y, hist_x, int(hist_y + dy), int(hist_x + dx)
  216. )
  217. draw.set_color(descs_img, (rows, cols), color, alpha=val)
  218. return descs, descs_img
  219. else:
  220. return descs