shape.py 7.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247
  1. import numbers
  2. import numpy as np
  3. from numpy.lib.stride_tricks import as_strided
  4. __all__ = ['view_as_blocks', 'view_as_windows']
  5. def view_as_blocks(arr_in, block_shape):
  6. """Block view of the input n-dimensional array (using re-striding).
  7. Blocks are non-overlapping views of the input array.
  8. Parameters
  9. ----------
  10. arr_in : ndarray, shape (M[, ...])
  11. Input array.
  12. block_shape : tuple
  13. The shape of the block. Each dimension must divide evenly into the
  14. corresponding dimensions of `arr_in`.
  15. Returns
  16. -------
  17. arr_out : ndarray
  18. Block view of the input array.
  19. Examples
  20. --------
  21. >>> import numpy as np
  22. >>> from skimage.util.shape import view_as_blocks
  23. >>> A = np.arange(4*4).reshape(4,4)
  24. >>> A
  25. array([[ 0, 1, 2, 3],
  26. [ 4, 5, 6, 7],
  27. [ 8, 9, 10, 11],
  28. [12, 13, 14, 15]])
  29. >>> B = view_as_blocks(A, block_shape=(2, 2))
  30. >>> B[0, 0]
  31. array([[0, 1],
  32. [4, 5]])
  33. >>> B[0, 1]
  34. array([[2, 3],
  35. [6, 7]])
  36. >>> B[1, 0, 1, 1]
  37. 13
  38. >>> A = np.arange(4*4*6).reshape(4,4,6)
  39. >>> A # doctest: +NORMALIZE_WHITESPACE
  40. array([[[ 0, 1, 2, 3, 4, 5],
  41. [ 6, 7, 8, 9, 10, 11],
  42. [12, 13, 14, 15, 16, 17],
  43. [18, 19, 20, 21, 22, 23]],
  44. [[24, 25, 26, 27, 28, 29],
  45. [30, 31, 32, 33, 34, 35],
  46. [36, 37, 38, 39, 40, 41],
  47. [42, 43, 44, 45, 46, 47]],
  48. [[48, 49, 50, 51, 52, 53],
  49. [54, 55, 56, 57, 58, 59],
  50. [60, 61, 62, 63, 64, 65],
  51. [66, 67, 68, 69, 70, 71]],
  52. [[72, 73, 74, 75, 76, 77],
  53. [78, 79, 80, 81, 82, 83],
  54. [84, 85, 86, 87, 88, 89],
  55. [90, 91, 92, 93, 94, 95]]])
  56. >>> B = view_as_blocks(A, block_shape=(1, 2, 2))
  57. >>> B.shape
  58. (4, 2, 3, 1, 2, 2)
  59. >>> B[2:, 0, 2] # doctest: +NORMALIZE_WHITESPACE
  60. array([[[[52, 53],
  61. [58, 59]]],
  62. [[[76, 77],
  63. [82, 83]]]])
  64. """
  65. if not isinstance(block_shape, tuple):
  66. raise TypeError('block needs to be a tuple')
  67. block_shape = np.array(block_shape)
  68. if (block_shape <= 0).any():
  69. raise ValueError("'block_shape' elements must be strictly positive")
  70. if block_shape.size != arr_in.ndim:
  71. raise ValueError("'block_shape' must have the same length " "as 'arr_in.shape'")
  72. arr_shape = np.array(arr_in.shape)
  73. if (arr_shape % block_shape).sum() != 0:
  74. raise ValueError("'block_shape' is not compatible with 'arr_in'")
  75. # -- restride the array to build the block view
  76. new_shape = tuple(arr_shape // block_shape) + tuple(block_shape)
  77. new_strides = tuple(arr_in.strides * block_shape) + arr_in.strides
  78. arr_out = as_strided(arr_in, shape=new_shape, strides=new_strides)
  79. return arr_out
  80. def view_as_windows(arr_in, window_shape, step=1):
  81. """Rolling window view of the input n-dimensional array.
  82. Windows are overlapping views of the input array, with adjacent windows
  83. shifted by a single row or column (or an index of a higher dimension).
  84. Parameters
  85. ----------
  86. arr_in : ndarray, shape (M[, ...])
  87. Input array.
  88. window_shape : integer or tuple of length arr_in.ndim
  89. Defines the shape of the elementary n-dimensional orthotope
  90. (better know as hyperrectangle [1]_) of the rolling window view.
  91. If an integer is given, the shape will be a hypercube of
  92. sidelength given by its value.
  93. step : integer or tuple of length arr_in.ndim
  94. Indicates step size at which extraction shall be performed.
  95. If integer is given, then the step is uniform in all dimensions.
  96. Returns
  97. -------
  98. arr_out : ndarray
  99. (rolling) window view of the input array.
  100. Notes
  101. -----
  102. One should be very careful with rolling views when it comes to
  103. memory usage. Indeed, although a 'view' has the same memory
  104. footprint as its base array, the actual array that emerges when this
  105. 'view' is used in a computation is generally a (much) larger array
  106. than the original, especially for 2-dimensional arrays and above.
  107. For example, let us consider a 3 dimensional array of size (100,
  108. 100, 100) of ``float64``. This array takes about 8*100**3 Bytes for
  109. storage which is just 8 MB. If one decides to build a rolling view
  110. on this array with a window of (3, 3, 3) the hypothetical size of
  111. the rolling view (if one was to reshape the view for example) would
  112. be 8*(100-3+1)**3*3**3 which is about 203 MB! The scaling becomes
  113. even worse as the dimension of the input array becomes larger.
  114. References
  115. ----------
  116. .. [1] https://en.wikipedia.org/wiki/Hyperrectangle
  117. Examples
  118. --------
  119. >>> import numpy as np
  120. >>> from skimage.util.shape import view_as_windows
  121. >>> A = np.arange(4*4).reshape(4,4)
  122. >>> A
  123. array([[ 0, 1, 2, 3],
  124. [ 4, 5, 6, 7],
  125. [ 8, 9, 10, 11],
  126. [12, 13, 14, 15]])
  127. >>> window_shape = (2, 2)
  128. >>> B = view_as_windows(A, window_shape)
  129. >>> B[0, 0]
  130. array([[0, 1],
  131. [4, 5]])
  132. >>> B[0, 1]
  133. array([[1, 2],
  134. [5, 6]])
  135. >>> A = np.arange(10)
  136. >>> A
  137. array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
  138. >>> window_shape = (3,)
  139. >>> B = view_as_windows(A, window_shape)
  140. >>> B.shape
  141. (8, 3)
  142. >>> B
  143. array([[0, 1, 2],
  144. [1, 2, 3],
  145. [2, 3, 4],
  146. [3, 4, 5],
  147. [4, 5, 6],
  148. [5, 6, 7],
  149. [6, 7, 8],
  150. [7, 8, 9]])
  151. >>> A = np.arange(5*4).reshape(5, 4)
  152. >>> A
  153. array([[ 0, 1, 2, 3],
  154. [ 4, 5, 6, 7],
  155. [ 8, 9, 10, 11],
  156. [12, 13, 14, 15],
  157. [16, 17, 18, 19]])
  158. >>> window_shape = (4, 3)
  159. >>> B = view_as_windows(A, window_shape)
  160. >>> B.shape
  161. (2, 2, 4, 3)
  162. >>> B # doctest: +NORMALIZE_WHITESPACE
  163. array([[[[ 0, 1, 2],
  164. [ 4, 5, 6],
  165. [ 8, 9, 10],
  166. [12, 13, 14]],
  167. [[ 1, 2, 3],
  168. [ 5, 6, 7],
  169. [ 9, 10, 11],
  170. [13, 14, 15]]],
  171. [[[ 4, 5, 6],
  172. [ 8, 9, 10],
  173. [12, 13, 14],
  174. [16, 17, 18]],
  175. [[ 5, 6, 7],
  176. [ 9, 10, 11],
  177. [13, 14, 15],
  178. [17, 18, 19]]]])
  179. """
  180. # -- basic checks on arguments
  181. if not isinstance(arr_in, np.ndarray):
  182. raise TypeError("`arr_in` must be a numpy ndarray")
  183. ndim = arr_in.ndim
  184. if isinstance(window_shape, numbers.Number):
  185. window_shape = (window_shape,) * ndim
  186. if not (len(window_shape) == ndim):
  187. raise ValueError("`window_shape` is incompatible with `arr_in.shape`")
  188. if isinstance(step, numbers.Number):
  189. if step < 1:
  190. raise ValueError("`step` must be >= 1")
  191. step = (step,) * ndim
  192. if len(step) != ndim:
  193. raise ValueError("`step` is incompatible with `arr_in.shape`")
  194. arr_shape = np.array(arr_in.shape)
  195. window_shape = np.array(window_shape, dtype=arr_shape.dtype)
  196. if ((arr_shape - window_shape) < 0).any():
  197. raise ValueError("`window_shape` is too large")
  198. if ((window_shape - 1) < 0).any():
  199. raise ValueError("`window_shape` is too small")
  200. # -- build rolling window view
  201. slices = tuple(slice(None, None, st) for st in step)
  202. window_strides = np.array(arr_in.strides)
  203. indexing_strides = arr_in[slices].strides
  204. win_indices_shape = (
  205. (np.array(arr_in.shape) - np.array(window_shape)) // np.array(step)
  206. ) + 1
  207. new_shape = tuple(list(win_indices_shape) + list(window_shape))
  208. strides = tuple(list(indexing_strides) + list(window_strides))
  209. arr_out = as_strided(arr_in, shape=new_shape, strides=strides)
  210. return arr_out