numpy-based filtergrid2d function (useful to apply filters on images in pure numpy)

filtergrid2d.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## scipy generic_filter like feature in pure numpy\n",
    "\n",
    "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",
    "\n",
    "For now it is only done for 2D arrays, but it should be possible to extend to nD arrays.\n",
    "\n",
    "The initial purpose was to do an image median filter which ignores NaNs, but this can be used for any kind of window-based filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "from numpy.lib import stride_tricks\n",
    "\n",
    "\n",
    "def filtergrid2d(data, window_size, mode='none', cval=0, flattenwindow=True):\n",
    "    \"\"\"Create an array of all filter windows of an image\n",
    "    \n",
    "    :param numpy.ndarray data: The array for which to create the grid\n",
    "    :param Union[int,List[int]] window_size:\n",
    "        Size of the window area for each dimensions or a\n",
    "        single int for window area with same size on each dimension\n",
    "    :param str mode: {'none', 'nearest', 'constant', 'reflect', 'mirror', 'wrap'}\n",
    "        Set the way to handle array borders.\n",
    "        Default is 'none' (output array shape is smaller than input array).\n",
    "    :param float cval:\n",
    "        Value used to fill borders when mode is 'constant'. Default is 0.\n",
    "    :param bool flattenwindow:\n",
    "        True (default) to return each window as 1D arrays\n",
    "        False to return windows keeping their shape\n",
    "    :return: Array of windows corresponding to each input array elements\n",
    "    :rtype: numpy.ndarray\n",
    "    \"\"\"\n",
    "    # Check input parameters\n",
    "    \n",
    "    if isinstance(window_size, int):\n",
    "        window_size = window_size, window_size\n",
    "    if window_size[0] <=0 or window_size[1] <= 0:\n",
    "        raise ValueError('Window size must be strictly positive')\n",
    "    \n",
    "    data = numpy.array(data, copy=False)\n",
    "    \n",
    "    if data.size == 0:\n",
    "        raise ValueError('zero-size array not supported')\n",
    "    \n",
    "    # Handle mode\n",
    "    \n",
    "    if mode == 'none':  # Use data as is provided\n",
    "        data = numpy.ascontiguousarray(data)\n",
    "        \n",
    "    else:  # Add borders to input data\n",
    "        input_data = data\n",
    "        \n",
    "        # Allocate an array with borders\n",
    "        data = numpy.empty((input_data.shape[0] + window_size[0] - 1,\n",
    "                            input_data.shape[1] + window_size[1] - 1),\n",
    "                           dtype=data.dtype)\n",
    "        start = window_size[0] // 2, window_size[1] // 2\n",
    "        end = start[0] + input_data.shape[0], start[1] + input_data.shape[1]\n",
    "        \n",
    "        data[start[0]:end[0], start[1]:end[1]] = input_data\n",
    "\n",
    "        if mode == 'constant':\n",
    "            data[:, :start[1]] = cval\n",
    "            data[:, end[1]:] = cval\n",
    "            data[:start[0], :] = cval\n",
    "            data[end[0]:, :] = cval\n",
    "            \n",
    "        elif mode == 'nearest':\n",
    "            # Fill borders\n",
    "            data[start[0]:end[0], :start[1]] = input_data[:, 0].reshape(-1, 1)\n",
    "            data[start[0]:end[0], end[1]:] = input_data[:, -1].reshape(-1, 1)\n",
    "            data[:start[0], start[1]:end[1]] = input_data[0, :]\n",
    "            data[end[0]:, start[1]:end[1]] = input_data[-1, :]\n",
    "            # Fill corners\n",
    "            data[:start[0], :start[1]] = input_data[0, 0]\n",
    "            data[:start[0], end[1]:] = input_data[0, -1]\n",
    "            data[end[0]:, :start[1]] = input_data[-1, 0]\n",
    "            data[end[0]:, end[1]:] = input_data[-1, -1]\n",
    "        \n",
    "        elif mode == 'reflect':\n",
    "            if start[0] > input_data.shape[0] or start[1] > input_data.shape[1]:\n",
    "                raise NotImplementedError(\n",
    "                    'window_size twice larger than data size not implemented')\n",
    "\n",
    "            end_offset = data.shape[0] - end[0], data.shape[1] - end[1]\n",
    "\n",
    "            # Fill borders\n",
    "            if start[1] > 0:\n",
    "                data[start[0]:end[0], :start[1]] = input_data[:, start[1]-1::-1]\n",
    "            data[start[0]:end[0], end[1]:] = input_data[:, -1:-1-end_offset[1]:-1]\n",
    "            if start[0] > 0:\n",
    "                data[:start[0], start[1]:end[1]] = input_data[start[0]-1::-1, :]\n",
    "            data[end[0]:, start[1]:end[1]] = input_data[-1:-1-end_offset[0]:-1, :]\n",
    "            # Fill corners\n",
    "            if start[0] > 0 and start[1] > 0:\n",
    "                data[:start[0], :start[1]] = input_data[start[0]-1::-1, start[1]-1::-1]\n",
    "            if start[0] > 0:\n",
    "                data[:start[0], end[1]:] = input_data[start[0]-1::-1, -1:-1-end_offset[1]:-1]\n",
    "            if start[1] > 0:\n",
    "                data[end[0]:, :start[1]] = input_data[-1:-1-end_offset[0]:-1, start[1]-1::-1]\n",
    "            data[end[0]:, end[1]:] = input_data[-1:-1-end_offset[0]:-1, -1:-1-end_offset[1]:-1]\n",
    "\n",
    "        elif mode == 'mirror':\n",
    "            raise NotImplementedError('mirror mode not implemented')\n",
    "        \n",
    "        elif mode == 'wrap':\n",
    "            raise NotImplementedError('wrap mode not implemented')\n",
    "            \n",
    "        else:\n",
    "            raise ValueError('Unsupported mode: %s' % str(mode))\n",
    "            \n",
    "    # Generate the grid\n",
    "    \n",
    "    itemsize = data.itemsize\n",
    "    dim0, dim1 = data.shape\n",
    "    \n",
    "    window_dim0 = min(window_size[0], dim0)\n",
    "    window_dim1 = min(window_size[1], dim1)\n",
    "\n",
    "    # Create a 4D view of the data that is convenient to apply the filter at once\n",
    "    # view uses memory overrides and overflows allocated memory...\n",
    "    view = stride_tricks.as_strided(\n",
    "        data,\n",
    "        shape=(dim0 + 1 - window_dim0, dim1, window_dim0, dim1),\n",
    "        strides=(itemsize * dim1, itemsize, itemsize * dim1, itemsize))\n",
    "        #numpy >=1.12 , writeable=False)\n",
    "\n",
    "    # Remove access to data that is not needed to apply filter\n",
    "    grid = view[:, :(dim1 + 1 - window_dim1), :, :window_dim1]\n",
    "\n",
    "    # Reshape array\n",
    "    shape = dim0 + 1 - window_dim0, dim1 + 1 - window_dim1\n",
    "    if flattenwindow:  # shape: (output dim0, output dim1, window size)\n",
    "        shape = shape + (window_dim0 * window_dim1,)\n",
    "    else:  # shape: (output dim0, output dim1, window dim0, window dim1)\n",
    "        shape = shape + (window_dim0, window_dim1)\n",
    "    grid = grid.reshape(*shape)\n",
    "    \n",
    "    return grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image:\n",
      "[[ 0  1  2  3]\n",
      " [ 4  5  6  7]\n",
      " [ 8  9 10 11]\n",
      " [12 13 14 15]]\n",
      "\n",
      "Filter grid:\n",
      "[[[[ 0  1  2]\n",
      "   [ 4  5  6]\n",
      "   [ 8  9 10]]\n",
      "\n",
      "  [[ 1  2  3]\n",
      "   [ 5  6  7]\n",
      "   [ 9 10 11]]]\n",
      "\n",
      "\n",
      " [[[ 4  5  6]\n",
      "   [ 8  9 10]\n",
      "   [12 13 14]]\n",
      "\n",
      "  [[ 5  6  7]\n",
      "   [ 9 10 11]\n",
      "   [13 14 15]]]]\n",
      "\n",
      "Filter grid (flatten window):\n",
      "[[[ 0  1  2  4  5  6  8  9 10]\n",
      "  [ 1  2  3  5  6  7  9 10 11]]\n",
      "\n",
      " [[ 4  5  6  8  9 10 12 13 14]\n",
      "  [ 5  6  7  9 10 11 13 14 15]]]\n"
     ]
    }
   ],
   "source": [
    "image = numpy.arange(16).reshape(4, 4)\n",
    "print('Image:')\n",
    "print(image)\n",
    "\n",
    "grid = filtergrid2d(image, window_size=(3, 3), flattenwindow=False)\n",
    "print('\\nFilter grid:')\n",
    "print(grid)\n",
    "\n",
    "flattengrid = filtergrid2d(image, window_size=(3, 3), flattenwindow=True)\n",
    "print('\\nFilter grid (flatten window):')\n",
    "print(flattengrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nanmedfilt2d(image, window_size=3):\n",
    "    return numpy.nanmedian(filtergrid2d(image, window_size), axis=-1)\n",
    "\n",
    "def medfilt2d(image, window_size=3):\n",
    "    return numpy.median(filtergrid2d(image, window_size), axis=-1)\n",
    "\n",
    "def maxfilt2d(image, window_size=3):\n",
    "    return numpy.max(filtergrid2d(image, window_size), axis=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basic tests: compare with generic_filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mode=\"nearest\" ...\n",
      "mode=\"constant\" ...\n",
      "mode=\"reflect\" ...\n",
      "    - window_size=18: NotImplementedError('window_size twice larger than data size not implemented',)\n"
     ]
    }
   ],
   "source": [
    "from scipy.ndimage.filters import generic_filter\n",
    "\n",
    "image = numpy.random.random(80).reshape(10, 8)\n",
    "function = numpy.mean\n",
    "\n",
    "for mode in ('nearest', 'constant', 'reflect'):\n",
    "    print('mode=\"%s\" ...' % mode)\n",
    "    for window_size in (1, 2, 3, 4, 5, (2, 3), (4, 5), (7, 5), 18):\n",
    "        ref = generic_filter(image, function, size=window_size, mode=mode)\n",
    "        try:\n",
    "            result = function(filtergrid2d(image, window_size, mode=mode), axis=-1)\n",
    "        except Exception as e:\n",
    "            print('    - window_size={}: {}'.format(window_size, repr(e)))\n",
    "        else:\n",
    "            if not numpy.all(numpy.equal(ref, result)):\n",
    "                print('    - window_size={}: {}'.format(window_size, 'Result mismatch'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Benchmark for median filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numpy: 1.15.1\n",
      "scipy: 1.1.0\n"
     ]
    }
   ],
   "source": [
    "import numpy\n",
    "import scipy\n",
    "\n",
    "print('numpy:', numpy.version.version)\n",
    "print('scipy:', scipy.version.version)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# benchmark images\n",
    "images = [numpy.random.random(500**2).reshape(500, 500), numpy.random.random(2048**2).reshape(2048, 2048)]\n",
    "\n",
    "# benchmark window size\n",
    "window_sizes = [(5, 5), (21, 21)]\n",
    "\n",
    "# benchmark conditions\n",
    "tests = [(image, window_size) for window_size in window_sizes for image in images]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "image shape (500, 500) window (5, 5) :\n",
      "181 ms ± 6.61 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "image shape (2048, 2048) window (5, 5) :\n",
      "2.87 s ± 39.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "image shape (500, 500) window (21, 21) :\n",
      "2.18 s ± 8.12 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "image shape (2048, 2048) window (21, 21) :\n",
      "36.7 s ± 346 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "# filtergrid2d + median, reflect mode\n",
    "\n",
    "for image, size in tests:\n",
    "    print('image shape', image.shape, 'window', size, ':')\n",
    "    %timeit numpy.median(filtergrid2d(image, size, mode='reflect'), axis=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "image shape (500, 500) window (5, 5) :\n",
      "Duration of a single run 12.62424898147583 s\n",
      "image shape (2048, 2048) window (5, 5) :\n",
      "Duration of a single run 205.4341971874237 s\n",
      "image shape (500, 500) window (21, 21) :\n",
      "Duration of a single run 14.224368810653687 s\n",
      "image shape (2048, 2048) window (21, 21) :\n",
      "Duration of a single run 235.9950499534607 s\n"
     ]
    }
   ],
   "source": [
    "# scipy generic filter + median, reflect mode\n",
    "# !!! slow !!!\n",
    "from scipy.ndimage.filters import generic_filter\n",
    "import time\n",
    "\n",
    "for image, size in tests:\n",
    "    print('image shape', image.shape, 'window', size, ':')\n",
    "    st = time.time()\n",
    "    result = generic_filter(image, numpy.median, size, mode='reflect')\n",
    "    print('Duration of a single run', time.time() - st, 's')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "image shape (500, 500) window (5, 5) :\n",
      "122 ms ± 5.54 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "image shape (2048, 2048) window (5, 5) :\n",
      "2.01 s ± 968 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "image shape (500, 500) window (21, 21) :\n",
      "1.68 s ± 193 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "image shape (2048, 2048) window (21, 21) :\n",
      "28.1 s ± 14.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "# scipy.ndimage.filters.median_filter, reflect mode\n",
    "from scipy.ndimage.filters import median_filter\n",
    "\n",
    "for image, size in tests:\n",
    "    print('image shape', image.shape, 'window', size, ':')\n",
    "    %timeit median_filter(image, size, mode='reflect')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}