pyramids.py 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408
  1. import math
  2. import numpy as np
  3. from .._shared.filters import gaussian
  4. from .._shared.utils import convert_to_float
  5. from ._warps import resize
  6. def _smooth(image, sigma, mode, cval, channel_axis):
  7. """Return image with each channel smoothed by the Gaussian filter."""
  8. smoothed = np.empty_like(image)
  9. # apply Gaussian filter to all channels independently
  10. if channel_axis is not None:
  11. # can rely on gaussian to insert a 0 entry at channel_axis
  12. channel_axis = channel_axis % image.ndim
  13. sigma = (sigma,) * (image.ndim - 1)
  14. else:
  15. channel_axis = None
  16. gaussian(
  17. image,
  18. sigma=sigma,
  19. out=smoothed,
  20. mode=mode,
  21. cval=cval,
  22. channel_axis=channel_axis,
  23. )
  24. return smoothed
  25. def _check_factor(factor):
  26. if factor <= 1:
  27. raise ValueError('scale factor must be greater than 1')
  28. def pyramid_reduce(
  29. image,
  30. downscale=2,
  31. sigma=None,
  32. order=1,
  33. mode='reflect',
  34. cval=0,
  35. preserve_range=False,
  36. *,
  37. channel_axis=None,
  38. ):
  39. """Smooth and then downsample image.
  40. Parameters
  41. ----------
  42. image : ndarray
  43. Input image.
  44. downscale : float, optional
  45. Downscale factor.
  46. sigma : float, optional
  47. Sigma for Gaussian filter. Default is `2 * downscale / 6.0` which
  48. corresponds to a filter mask twice the size of the scale factor that
  49. covers more than 99% of the Gaussian distribution.
  50. order : int, optional
  51. Order of splines used in interpolation of downsampling. See
  52. `skimage.transform.warp` for detail.
  53. mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
  54. The mode parameter determines how the array borders are handled, where
  55. cval is the value when mode is equal to 'constant'.
  56. cval : float, optional
  57. Value to fill past edges of input if mode is 'constant'.
  58. preserve_range : bool, optional
  59. Whether to keep the original range of values. Otherwise, the input
  60. image is converted according to the conventions of `img_as_float`.
  61. Also see https://scikit-image.org/docs/dev/user_guide/data_types.html
  62. channel_axis : int or None, optional
  63. If None, the image is assumed to be a grayscale (single channel) image.
  64. Otherwise, this parameter indicates which axis of the array corresponds
  65. to channels.
  66. .. versionadded:: 0.19
  67. ``channel_axis`` was added in 0.19.
  68. Returns
  69. -------
  70. out : array
  71. Smoothed and downsampled float image.
  72. References
  73. ----------
  74. .. [1] http://persci.mit.edu/pub_pdfs/pyramid83.pdf
  75. """
  76. _check_factor(downscale)
  77. image = convert_to_float(image, preserve_range)
  78. if channel_axis is not None:
  79. channel_axis = channel_axis % image.ndim
  80. out_shape = tuple(
  81. math.ceil(d / float(downscale)) if ax != channel_axis else d
  82. for ax, d in enumerate(image.shape)
  83. )
  84. else:
  85. out_shape = tuple(math.ceil(d / float(downscale)) for d in image.shape)
  86. if sigma is None:
  87. # automatically determine sigma which covers > 99% of distribution
  88. sigma = 2 * downscale / 6.0
  89. smoothed = _smooth(image, sigma, mode, cval, channel_axis)
  90. out = resize(
  91. smoothed, out_shape, order=order, mode=mode, cval=cval, anti_aliasing=False
  92. )
  93. return out
  94. def pyramid_expand(
  95. image,
  96. upscale=2,
  97. sigma=None,
  98. order=1,
  99. mode='reflect',
  100. cval=0,
  101. preserve_range=False,
  102. *,
  103. channel_axis=None,
  104. ):
  105. """Upsample and then smooth image.
  106. Parameters
  107. ----------
  108. image : ndarray
  109. Input image.
  110. upscale : float, optional
  111. Upscale factor.
  112. sigma : float, optional
  113. Sigma for Gaussian filter. Default is `2 * upscale / 6.0` which
  114. corresponds to a filter mask twice the size of the scale factor that
  115. covers more than 99% of the Gaussian distribution.
  116. order : int, optional
  117. Order of splines used in interpolation of upsampling. See
  118. `skimage.transform.warp` for detail.
  119. mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
  120. The mode parameter determines how the array borders are handled, where
  121. cval is the value when mode is equal to 'constant'.
  122. cval : float, optional
  123. Value to fill past edges of input if mode is 'constant'.
  124. preserve_range : bool, optional
  125. Whether to keep the original range of values. Otherwise, the input
  126. image is converted according to the conventions of `img_as_float`.
  127. Also see https://scikit-image.org/docs/dev/user_guide/data_types.html
  128. channel_axis : int or None, optional
  129. If None, the image is assumed to be a grayscale (single channel) image.
  130. Otherwise, this parameter indicates which axis of the array corresponds
  131. to channels.
  132. .. versionadded:: 0.19
  133. ``channel_axis`` was added in 0.19.
  134. Returns
  135. -------
  136. out : array
  137. Upsampled and smoothed float image.
  138. References
  139. ----------
  140. .. [1] http://persci.mit.edu/pub_pdfs/pyramid83.pdf
  141. """
  142. _check_factor(upscale)
  143. image = convert_to_float(image, preserve_range)
  144. if channel_axis is not None:
  145. channel_axis = channel_axis % image.ndim
  146. out_shape = tuple(
  147. math.ceil(upscale * d) if ax != channel_axis else d
  148. for ax, d in enumerate(image.shape)
  149. )
  150. else:
  151. out_shape = tuple(math.ceil(upscale * d) for d in image.shape)
  152. if sigma is None:
  153. # automatically determine sigma which covers > 99% of distribution
  154. sigma = 2 * upscale / 6.0
  155. resized = resize(
  156. image, out_shape, order=order, mode=mode, cval=cval, anti_aliasing=False
  157. )
  158. out = _smooth(resized, sigma, mode, cval, channel_axis)
  159. return out
  160. def pyramid_gaussian(
  161. image,
  162. max_layer=-1,
  163. downscale=2,
  164. sigma=None,
  165. order=1,
  166. mode='reflect',
  167. cval=0,
  168. preserve_range=False,
  169. *,
  170. channel_axis=None,
  171. ):
  172. """Yield images of the Gaussian pyramid formed by the input image.
  173. Recursively applies the `pyramid_reduce` function to the image, and yields
  174. the downscaled images.
  175. Note that the first image of the pyramid will be the original, unscaled
  176. image. The total number of images is `max_layer + 1`. In case all layers
  177. are computed, the last image is either a one-pixel image or the image where
  178. the reduction does not change its shape.
  179. Parameters
  180. ----------
  181. image : ndarray
  182. Input image.
  183. max_layer : int, optional
  184. Number of layers for the pyramid. 0th layer is the original image.
  185. Default is -1 which builds all possible layers.
  186. downscale : float, optional
  187. Downscale factor.
  188. sigma : float, optional
  189. Sigma for Gaussian filter. Default is `2 * downscale / 6.0` which
  190. corresponds to a filter mask twice the size of the scale factor that
  191. covers more than 99% of the Gaussian distribution.
  192. order : int, optional
  193. Order of splines used in interpolation of downsampling. See
  194. `skimage.transform.warp` for detail.
  195. mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
  196. The mode parameter determines how the array borders are handled, where
  197. cval is the value when mode is equal to 'constant'.
  198. cval : float, optional
  199. Value to fill past edges of input if mode is 'constant'.
  200. preserve_range : bool, optional
  201. Whether to keep the original range of values. Otherwise, the input
  202. image is converted according to the conventions of `img_as_float`.
  203. Also see https://scikit-image.org/docs/dev/user_guide/data_types.html
  204. channel_axis : int or None, optional
  205. If None, the image is assumed to be a grayscale (single channel) image.
  206. Otherwise, this parameter indicates which axis of the array corresponds
  207. to channels.
  208. .. versionadded:: 0.19
  209. ``channel_axis`` was added in 0.19.
  210. Returns
  211. -------
  212. pyramid : generator
  213. Generator yielding pyramid layers as float images.
  214. References
  215. ----------
  216. .. [1] http://persci.mit.edu/pub_pdfs/pyramid83.pdf
  217. """
  218. _check_factor(downscale)
  219. # cast to float for consistent data type in pyramid
  220. image = convert_to_float(image, preserve_range)
  221. layer = 0
  222. current_shape = image.shape
  223. prev_layer_image = image
  224. yield image
  225. # build downsampled images until max_layer is reached or downscale process
  226. # does not change image size
  227. while layer != max_layer:
  228. layer += 1
  229. layer_image = pyramid_reduce(
  230. prev_layer_image,
  231. downscale,
  232. sigma,
  233. order,
  234. mode,
  235. cval,
  236. channel_axis=channel_axis,
  237. )
  238. prev_shape = current_shape
  239. prev_layer_image = layer_image
  240. current_shape = layer_image.shape
  241. # no change to previous pyramid layer
  242. if current_shape == prev_shape:
  243. break
  244. yield layer_image
  245. def pyramid_laplacian(
  246. image,
  247. max_layer=-1,
  248. downscale=2,
  249. sigma=None,
  250. order=1,
  251. mode='reflect',
  252. cval=0,
  253. preserve_range=False,
  254. *,
  255. channel_axis=None,
  256. ):
  257. """Yield images of the laplacian pyramid formed by the input image.
  258. Each layer contains the difference between the downsampled and the
  259. downsampled, smoothed image::
  260. layer = resize(prev_layer) - smooth(resize(prev_layer))
  261. Note that the first image of the pyramid will be the difference between the
  262. original, unscaled image and its smoothed version. The total number of
  263. images is `max_layer + 1`. In case all layers are computed, the last image
  264. is either a one-pixel image or the image where the reduction does not
  265. change its shape.
  266. Parameters
  267. ----------
  268. image : ndarray
  269. Input image.
  270. max_layer : int, optional
  271. Number of layers for the pyramid. 0th layer is the original image.
  272. Default is -1 which builds all possible layers.
  273. downscale : float, optional
  274. Downscale factor.
  275. sigma : float, optional
  276. Sigma for Gaussian filter. Default is `2 * downscale / 6.0` which
  277. corresponds to a filter mask twice the size of the scale factor that
  278. covers more than 99% of the Gaussian distribution.
  279. order : int, optional
  280. Order of splines used in interpolation of downsampling. See
  281. `skimage.transform.warp` for detail.
  282. mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
  283. The mode parameter determines how the array borders are handled, where
  284. cval is the value when mode is equal to 'constant'.
  285. cval : float, optional
  286. Value to fill past edges of input if mode is 'constant'.
  287. preserve_range : bool, optional
  288. Whether to keep the original range of values. Otherwise, the input
  289. image is converted according to the conventions of `img_as_float`.
  290. Also see https://scikit-image.org/docs/dev/user_guide/data_types.html
  291. channel_axis : int or None, optional
  292. If None, the image is assumed to be a grayscale (single channel) image.
  293. Otherwise, this parameter indicates which axis of the array corresponds
  294. to channels.
  295. .. versionadded:: 0.19
  296. ``channel_axis`` was added in 0.19.
  297. Returns
  298. -------
  299. pyramid : generator
  300. Generator yielding pyramid layers as float images.
  301. References
  302. ----------
  303. .. [1] http://persci.mit.edu/pub_pdfs/pyramid83.pdf
  304. .. [2] http://sepwww.stanford.edu/data/media/public/sep/morgan/texturematch/paper_html/node3.html
  305. """
  306. _check_factor(downscale)
  307. # cast to float for consistent data type in pyramid
  308. image = convert_to_float(image, preserve_range)
  309. if sigma is None:
  310. # automatically determine sigma which covers > 99% of distribution
  311. sigma = 2 * downscale / 6.0
  312. current_shape = image.shape
  313. smoothed_image = _smooth(image, sigma, mode, cval, channel_axis)
  314. yield image - smoothed_image
  315. if channel_axis is not None:
  316. channel_axis = channel_axis % image.ndim
  317. shape_without_channels = list(current_shape)
  318. shape_without_channels.pop(channel_axis)
  319. shape_without_channels = tuple(shape_without_channels)
  320. else:
  321. shape_without_channels = current_shape
  322. # build downsampled images until max_layer is reached or downscale process
  323. # does not change image size
  324. if max_layer == -1:
  325. max_layer = math.ceil(math.log(max(shape_without_channels), downscale))
  326. for layer in range(max_layer):
  327. if channel_axis is not None:
  328. out_shape = tuple(
  329. math.ceil(d / float(downscale)) if ax != channel_axis else d
  330. for ax, d in enumerate(current_shape)
  331. )
  332. else:
  333. out_shape = tuple(math.ceil(d / float(downscale)) for d in current_shape)
  334. resized_image = resize(
  335. smoothed_image,
  336. out_shape,
  337. order=order,
  338. mode=mode,
  339. cval=cval,
  340. anti_aliasing=False,
  341. )
  342. smoothed_image = _smooth(resized_image, sigma, mode, cval, channel_axis)
  343. current_shape = resized_image.shape
  344. yield resized_image - smoothed_image