fourier.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306
  1. # Copyright (C) 2003-2005 Peter J. Verveer
  2. #
  3. # Redistribution and use in source and binary forms, with or without
  4. # modification, are permitted provided that the following conditions
  5. # are met:
  6. #
  7. # 1. Redistributions of source code must retain the above copyright
  8. # notice, this list of conditions and the following disclaimer.
  9. #
  10. # 2. Redistributions in binary form must reproduce the above
  11. # copyright notice, this list of conditions and the following
  12. # disclaimer in the documentation and/or other materials provided
  13. # with the distribution.
  14. #
  15. # 3. The name of the author may not be used to endorse or promote
  16. # products derived from this software without specific prior
  17. # written permission.
  18. #
  19. # THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
  20. # OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  21. # WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  22. # ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
  23. # DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  24. # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
  25. # GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  26. # INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
  27. # WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  28. # NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  29. # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  30. from __future__ import division, print_function, absolute_import
  31. import numpy
  32. from . import _ni_support
  33. from . import _nd_image
  34. __all__ = ['fourier_gaussian', 'fourier_uniform', 'fourier_ellipsoid',
  35. 'fourier_shift']
  36. def _get_output_fourier(output, input):
  37. if output is None:
  38. if input.dtype.type in [numpy.complex64, numpy.complex128,
  39. numpy.float32]:
  40. output = numpy.zeros(input.shape, dtype=input.dtype)
  41. else:
  42. output = numpy.zeros(input.shape, dtype=numpy.float64)
  43. elif type(output) is type:
  44. if output not in [numpy.complex64, numpy.complex128,
  45. numpy.float32, numpy.float64]:
  46. raise RuntimeError("output type not supported")
  47. output = numpy.zeros(input.shape, dtype=output)
  48. elif output.shape != input.shape:
  49. raise RuntimeError("output shape not correct")
  50. return output
  51. def _get_output_fourier_complex(output, input):
  52. if output is None:
  53. if input.dtype.type in [numpy.complex64, numpy.complex128]:
  54. output = numpy.zeros(input.shape, dtype=input.dtype)
  55. else:
  56. output = numpy.zeros(input.shape, dtype=numpy.complex128)
  57. elif type(output) is type:
  58. if output not in [numpy.complex64, numpy.complex128]:
  59. raise RuntimeError("output type not supported")
  60. output = numpy.zeros(input.shape, dtype=output)
  61. elif output.shape != input.shape:
  62. raise RuntimeError("output shape not correct")
  63. return output
  64. def fourier_gaussian(input, sigma, n=-1, axis=-1, output=None):
  65. """
  66. Multi-dimensional Gaussian fourier filter.
  67. The array is multiplied with the fourier transform of a Gaussian
  68. kernel.
  69. Parameters
  70. ----------
  71. input : array_like
  72. The input array.
  73. sigma : float or sequence
  74. The sigma of the Gaussian kernel. If a float, `sigma` is the same for
  75. all axes. If a sequence, `sigma` has to contain one value for each
  76. axis.
  77. n : int, optional
  78. If `n` is negative (default), then the input is assumed to be the
  79. result of a complex fft.
  80. If `n` is larger than or equal to zero, the input is assumed to be the
  81. result of a real fft, and `n` gives the length of the array before
  82. transformation along the real transform direction.
  83. axis : int, optional
  84. The axis of the real transform.
  85. output : ndarray, optional
  86. If given, the result of filtering the input is placed in this array.
  87. None is returned in this case.
  88. Returns
  89. -------
  90. fourier_gaussian : ndarray
  91. The filtered input.
  92. Examples
  93. --------
  94. >>> from scipy import ndimage, misc
  95. >>> import numpy.fft
  96. >>> import matplotlib.pyplot as plt
  97. >>> fig, (ax1, ax2) = plt.subplots(1, 2)
  98. >>> plt.gray() # show the filtered result in grayscale
  99. >>> ascent = misc.ascent()
  100. >>> input_ = numpy.fft.fft2(ascent)
  101. >>> result = ndimage.fourier_gaussian(input_, sigma=4)
  102. >>> result = numpy.fft.ifft2(result)
  103. >>> ax1.imshow(ascent)
  104. >>> ax2.imshow(result.real) # the imaginary part is an artifact
  105. >>> plt.show()
  106. """
  107. input = numpy.asarray(input)
  108. output = _get_output_fourier(output, input)
  109. axis = _ni_support._check_axis(axis, input.ndim)
  110. sigmas = _ni_support._normalize_sequence(sigma, input.ndim)
  111. sigmas = numpy.asarray(sigmas, dtype=numpy.float64)
  112. if not sigmas.flags.contiguous:
  113. sigmas = sigmas.copy()
  114. _nd_image.fourier_filter(input, sigmas, n, axis, output, 0)
  115. return output
  116. def fourier_uniform(input, size, n=-1, axis=-1, output=None):
  117. """
  118. Multi-dimensional uniform fourier filter.
  119. The array is multiplied with the fourier transform of a box of given
  120. size.
  121. Parameters
  122. ----------
  123. input : array_like
  124. The input array.
  125. size : float or sequence
  126. The size of the box used for filtering.
  127. If a float, `size` is the same for all axes. If a sequence, `size` has
  128. to contain one value for each axis.
  129. n : int, optional
  130. If `n` is negative (default), then the input is assumed to be the
  131. result of a complex fft.
  132. If `n` is larger than or equal to zero, the input is assumed to be the
  133. result of a real fft, and `n` gives the length of the array before
  134. transformation along the real transform direction.
  135. axis : int, optional
  136. The axis of the real transform.
  137. output : ndarray, optional
  138. If given, the result of filtering the input is placed in this array.
  139. None is returned in this case.
  140. Returns
  141. -------
  142. fourier_uniform : ndarray
  143. The filtered input.
  144. Examples
  145. --------
  146. >>> from scipy import ndimage, misc
  147. >>> import numpy.fft
  148. >>> import matplotlib.pyplot as plt
  149. >>> fig, (ax1, ax2) = plt.subplots(1, 2)
  150. >>> plt.gray() # show the filtered result in grayscale
  151. >>> ascent = misc.ascent()
  152. >>> input_ = numpy.fft.fft2(ascent)
  153. >>> result = ndimage.fourier_uniform(input_, size=20)
  154. >>> result = numpy.fft.ifft2(result)
  155. >>> ax1.imshow(ascent)
  156. >>> ax2.imshow(result.real) # the imaginary part is an artifact
  157. >>> plt.show()
  158. """
  159. input = numpy.asarray(input)
  160. output = _get_output_fourier(output, input)
  161. axis = _ni_support._check_axis(axis, input.ndim)
  162. sizes = _ni_support._normalize_sequence(size, input.ndim)
  163. sizes = numpy.asarray(sizes, dtype=numpy.float64)
  164. if not sizes.flags.contiguous:
  165. sizes = sizes.copy()
  166. _nd_image.fourier_filter(input, sizes, n, axis, output, 1)
  167. return output
  168. def fourier_ellipsoid(input, size, n=-1, axis=-1, output=None):
  169. """
  170. Multi-dimensional ellipsoid fourier filter.
  171. The array is multiplied with the fourier transform of a ellipsoid of
  172. given sizes.
  173. Parameters
  174. ----------
  175. input : array_like
  176. The input array.
  177. size : float or sequence
  178. The size of the box used for filtering.
  179. If a float, `size` is the same for all axes. If a sequence, `size` has
  180. to contain one value for each axis.
  181. n : int, optional
  182. If `n` is negative (default), then the input is assumed to be the
  183. result of a complex fft.
  184. If `n` is larger than or equal to zero, the input is assumed to be the
  185. result of a real fft, and `n` gives the length of the array before
  186. transformation along the real transform direction.
  187. axis : int, optional
  188. The axis of the real transform.
  189. output : ndarray, optional
  190. If given, the result of filtering the input is placed in this array.
  191. None is returned in this case.
  192. Returns
  193. -------
  194. fourier_ellipsoid : ndarray
  195. The filtered input.
  196. Notes
  197. -----
  198. This function is implemented for arrays of rank 1, 2, or 3.
  199. Examples
  200. --------
  201. >>> from scipy import ndimage, misc
  202. >>> import numpy.fft
  203. >>> import matplotlib.pyplot as plt
  204. >>> fig, (ax1, ax2) = plt.subplots(1, 2)
  205. >>> plt.gray() # show the filtered result in grayscale
  206. >>> ascent = misc.ascent()
  207. >>> input_ = numpy.fft.fft2(ascent)
  208. >>> result = ndimage.fourier_ellipsoid(input_, size=20)
  209. >>> result = numpy.fft.ifft2(result)
  210. >>> ax1.imshow(ascent)
  211. >>> ax2.imshow(result.real) # the imaginary part is an artifact
  212. >>> plt.show()
  213. """
  214. input = numpy.asarray(input)
  215. output = _get_output_fourier(output, input)
  216. axis = _ni_support._check_axis(axis, input.ndim)
  217. sizes = _ni_support._normalize_sequence(size, input.ndim)
  218. sizes = numpy.asarray(sizes, dtype=numpy.float64)
  219. if not sizes.flags.contiguous:
  220. sizes = sizes.copy()
  221. _nd_image.fourier_filter(input, sizes, n, axis, output, 2)
  222. return output
  223. def fourier_shift(input, shift, n=-1, axis=-1, output=None):
  224. """
  225. Multi-dimensional fourier shift filter.
  226. The array is multiplied with the fourier transform of a shift operation.
  227. Parameters
  228. ----------
  229. input : array_like
  230. The input array.
  231. shift : float or sequence
  232. The size of the box used for filtering.
  233. If a float, `shift` is the same for all axes. If a sequence, `shift`
  234. has to contain one value for each axis.
  235. n : int, optional
  236. If `n` is negative (default), then the input is assumed to be the
  237. result of a complex fft.
  238. If `n` is larger than or equal to zero, the input is assumed to be the
  239. result of a real fft, and `n` gives the length of the array before
  240. transformation along the real transform direction.
  241. axis : int, optional
  242. The axis of the real transform.
  243. output : ndarray, optional
  244. If given, the result of shifting the input is placed in this array.
  245. None is returned in this case.
  246. Returns
  247. -------
  248. fourier_shift : ndarray
  249. The shifted input.
  250. Examples
  251. --------
  252. >>> from scipy import ndimage, misc
  253. >>> import matplotlib.pyplot as plt
  254. >>> import numpy.fft
  255. >>> fig, (ax1, ax2) = plt.subplots(1, 2)
  256. >>> plt.gray() # show the filtered result in grayscale
  257. >>> ascent = misc.ascent()
  258. >>> input_ = numpy.fft.fft2(ascent)
  259. >>> result = ndimage.fourier_shift(input_, shift=200)
  260. >>> result = numpy.fft.ifft2(result)
  261. >>> ax1.imshow(ascent)
  262. >>> ax2.imshow(result.real) # the imaginary part is an artifact
  263. >>> plt.show()
  264. """
  265. input = numpy.asarray(input)
  266. output = _get_output_fourier_complex(output, input)
  267. axis = _ni_support._check_axis(axis, input.ndim)
  268. shifts = _ni_support._normalize_sequence(shift, input.ndim)
  269. shifts = numpy.asarray(shifts, dtype=numpy.float64)
  270. if not shifts.flags.contiguous:
  271. shifts = shifts.copy()
  272. _nd_image.fourier_shift(input, shifts, n, axis, output)
  273. return output