API tools faq
paste
Guest User

Untitled

a guest
Jul 16th, 2018
77
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 8.27 KB | None | 0 0
raw download clone embed print report
  1. {
  2. "cells": [
  3. {
  4. "cell_type": "markdown",
  5. "metadata": {},
  6. "source": [
  7. "## scipy generic_filter like feature in pure numpy\n",
  8. "\n",
  9. "The idea is to create a single array from an image that provides the elements of a window for each pixel in the input image (function filtergrid2d).\n",
  10. "\n",
  11. "For now it is only done for 2D arrays, but it should be possible to extend to nD arrays.\n",
  12. "\n",
  13. "The initial purpose was to do an image median filter which ignores NaNs, but this can be used for any kind of filter based on window."
  14. ]
  15. },
  16. {
  17. "cell_type": "code",
  18. "execution_count": 19,
  19. "metadata": {},
  20. "outputs": [],
  21. "source": [
  22. "import numpy\n",
  23. "from numpy.lib import stride_tricks\n",
  24. "\n",
  25. "\n",
  26. "def filtergrid2d(data, size):\n",
  27. " \"\"\"Create a 3D array of all filter windows of an image\n",
  28. " \n",
  29. " :param numpy.ndarray data: The array for which to create the grid\n",
  30. " :param Union[int,List[int]] size:\n",
  31. " Size of the window area for each dimensions or a\n",
  32. " single int for window area with same size on each dimension\n",
  33. " :return: Array of windows corresponding to each input array elements\n",
  34. " (initial dimensions, window_size)\n",
  35. " :rtype: numpy.ndarray\n",
  36. " \"\"\"\n",
  37. " data = numpy.ascontiguousarray(data)\n",
  38. " itemsize = data.itemsize\n",
  39. " height, width = data.shape\n",
  40. "\n",
  41. " if isinstance(size, int):\n",
  42. " size = size, size\n",
  43. " kernel_height = min(size[0], height)\n",
  44. " kernel_width = min(size[1], width)\n",
  45. "\n",
  46. " # Create a 4D view of the data that is convenient to apply the filter at once\n",
  47. " # view uses memory overrides and overflows allocated memory...\n",
  48. " view = stride_tricks.as_strided(\n",
  49. " data,\n",
  50. " shape=(height + 1 - kernel_height, width, kernel_height, width),\n",
  51. " strides=(itemsize * width, itemsize, itemsize * width, itemsize))\n",
  52. " #numpy >=1.12 , writeable=False)\n",
  53. "\n",
  54. " # Remove access to data that is not needed to apply filter\n",
  55. " grid = view[:, :(width + 1 - kernel_width), :, :kernel_width]\n",
  56. " # Reshape array to (output height, output width, kernel size)\n",
  57. " grid = grid.reshape(height + 1 - kernel_height,\n",
  58. " width + 1 - kernel_width,\n",
  59. " kernel_height * kernel_width) # alternative kernel_height, kernel_width\n",
  60. " return grid"
  61. ]
  62. },
  63. {
  64. "cell_type": "markdown",
  65. "metadata": {},
  66. "source": [
  67. "### Example of use"
  68. ]
  69. },
  70. {
  71. "cell_type": "code",
  72. "execution_count": 20,
  73. "metadata": {},
  74. "outputs": [],
  75. "source": [
  76. "def nanmedfilt2d(image, size=(3, 3)):\n",
  77. " return numpy.nanmedian(filtergrid2d(image, size), axis=-1)\n",
  78. "\n",
  79. "def medfilt2d(image, size=(3, 3)):\n",
  80. " return numpy.median(filtergrid2d(image, size), axis=-1)"
  81. ]
  82. },
  83. {
  84. "cell_type": "markdown",
  85. "metadata": {},
  86. "source": [
  87. "### Benchmark"
  88. ]
  89. },
  90. {
  91. "cell_type": "code",
  92. "execution_count": 27,
  93. "metadata": {},
  94. "outputs": [],
  95. "source": [
  96. "# benchmark images\n",
  97. "images = [numpy.random.random(500**2).reshape(500, 500), numpy.random.random(2048**2).reshape(2048, 2048)]\n",
  98. "\n",
  99. "# benchmark window size\n",
  100. "sizes = [(5, 5), (21, 21)]\n",
  101. "\n",
  102. "# benchmark conditions\n",
  103. "tests = [(image, size) for size in sizes for image in images]"
  104. ]
  105. },
  106. {
  107. "cell_type": "code",
  108. "execution_count": 28,
  109. "metadata": {},
  110. "outputs": [
  111. {
  112. "name": "stdout",
  113. "output_type": "stream",
  114. "text": [
  115. "image shape (500, 500) window (5, 5) :\n",
  116. "120 ms ± 1.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
  117. "image shape (2048, 2048) window (5, 5) :\n",
  118. "1.98 s ± 9.85 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
  119. "image shape (500, 500) window (21, 21) :\n",
  120. "1.48 s ± 863 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
  121. "image shape (2048, 2048) window (21, 21) :\n",
  122. "26.6 s ± 76.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
  123. ]
  124. }
  125. ],
  126. "source": [
  127. "# filtergrid2d\n",
  128. "\n",
  129. "from scipy.signal import medfilt2d as scipy_medfilt2d\n",
  130. "\n",
  131. "for image, size in tests:\n",
  132. " print('image shape', image.shape, 'window', size, ':')\n",
  133. " # Do not compare borders\n",
  134. " assert(numpy.abs(scipy_medfilt2d(image, size)[size[0]//2:-(size[0]//2), size[1]//2:-(size[1]//2)] -\n",
  135. " medfilt2d(image, size)).max() == 0.)\n",
  136. "\n",
  137. " %timeit numpy.median(filtergrid2d(image, size), axis=-1)"
  138. ]
  139. },
  140. {
  141. "cell_type": "code",
  142. "execution_count": 29,
  143. "metadata": {},
  144. "outputs": [
  145. {
  146. "name": "stdout",
  147. "output_type": "stream",
  148. "text": [
  149. "image shape (500, 500) window (5, 5) :\n",
  150. "75.1 ms ± 434 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
  151. "image shape (2048, 2048) window (5, 5) :\n",
  152. "1.27 s ± 6.88 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
  153. "image shape (500, 500) window (21, 21) :\n",
  154. "1.08 s ± 9.87 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
  155. "image shape (2048, 2048) window (21, 21) :\n",
  156. "18.5 s ± 223 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
  157. ]
  158. }
  159. ],
  160. "source": [
  161. "# scipy.signal.medfilt2d benchmark\n",
  162. "from scipy.signal import medfilt2d as scipy_medfilt2d\n",
  163. "\n",
  164. "for image, size in tests:\n",
  165. " print('image shape', image.shape, 'window', size, ':')\n",
  166. " %timeit scipy_medfilt2d(image, size)"
  167. ]
  168. },
  169. {
  170. "cell_type": "code",
  171. "execution_count": 30,
  172. "metadata": {},
  173. "outputs": [
  174. {
  175. "name": "stdout",
  176. "output_type": "stream",
  177. "text": [
  178. "image shape (500, 500) window (5, 5) :\n",
  179. "85.4 ms ± 202 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
  180. "image shape (2048, 2048) window (5, 5) :\n",
  181. "1.43 s ± 1.46 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
  182. "image shape (500, 500) window (21, 21) :\n",
  183. "1.11 s ± 897 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
  184. "image shape (2048, 2048) window (21, 21) :\n",
  185. "19 s ± 95.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
  186. ]
  187. }
  188. ],
  189. "source": [
  190. "# scipy.ndimage.filters.median_filter\n",
  191. "from scipy.ndimage.filters import median_filter\n",
  192. "\n",
  193. "for image, size in tests:\n",
  194. " print('image shape', image.shape, 'window', size, ':')\n",
  195. " %timeit median_filter(image, size, mode='nearest')"
  196. ]
  197. },
  198. {
  199. "cell_type": "code",
  200. "execution_count": 31,
  201. "metadata": {},
  202. "outputs": [
  203. {
  204. "name": "stdout",
  205. "output_type": "stream",
  206. "text": [
  207. "image shape (500, 500) window (5, 5) :\n",
  208. "Duration of a single run 10.746914625167847 s\n",
  209. "image shape (2048, 2048) window (5, 5) :\n",
  210. "Duration of a single run 173.32792830467224 s\n",
  211. "image shape (500, 500) window (21, 21) :\n",
  212. "Duration of a single run 11.813754081726074 s\n",
  213. "image shape (2048, 2048) window (21, 21) :\n",
  214. "Duration of a single run 199.75266814231873 s\n"
  215. ]
  216. }
  217. ],
  218. "source": [
  219. "# scipy generic filter + nanmedian\n",
  220. "# !!! slow !!!\n",
  221. "from scipy.ndimage.filters import generic_filter\n",
  222. "import time\n",
  223. "\n",
  224. "for image, size in tests:\n",
  225. " print('image shape', image.shape, 'window', size, ':')\n",
  226. " st = time.time()\n",
  227. " generic_filter(image, numpy.nanmedian, size)\n",
  228. " print('Duration of a single run', time.time() - st, 's')"
  229. ]
  230. },
  231. {
  232. "cell_type": "code",
  233. "execution_count": null,
  234. "metadata": {},
  235. "outputs": [],
  236. "source": []
  237. }
  238. ],
  239. "metadata": {
  240. "kernelspec": {
  241. "display_name": "Python 3",
  242. "language": "python",
  243. "name": "python3"
  244. },
  245. "language_info": {
  246. "codemirror_mode": {
  247. "name": "ipython",
  248. "version": 3
  249. },
  250. "file_extension": ".py",
  251. "mimetype": "text/x-python",
  252. "name": "python",
  253. "nbconvert_exporter": "python",
  254. "pygments_lexer": "ipython3",
  255. "version": "3.4.2"
  256. }
  257. },
  258. "nbformat": 4,
  259. "nbformat_minor": 2
  260. }
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!