{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "6223a8bf",
            "cell_type": "markdown",
            "source": "<div class='alert alert-warning'>\n\nSciPy's interactive examples with Jupyterlite are experimental and may not always work as expected. Execution of cells containing imports may result in large downloads (up to 60MB of content for the first import from SciPy). Load times when importing from SciPy may take roughly 10-20 seconds. If you notice any problems, feel free to open an [issue](https://github.com/scipy/scipy/issues/new/choose).\n\n</div>",
            "metadata": {}
        },
        {
            "id": "2b5d9e56",
            "cell_type": "markdown",
            "source": "Suppose we wish to perform a median filter with even window size on a ``float16``\nimage. Furthermore, the image has NaNs that we wish to be ignored (and effectively\nremoved by the filter). `median_filter` does not support ``float16`` data, its\nbehavior when NaNs are present is not defined, and for even window sizes, it does\nnot return the usual sample median - the average of the two middle elements. This\nwould be an excellent use case for `vectorized_filter` with\n``function=np.nanmedian``, which supports the required interface: it accepts a\ndata array of any shape as the first positional argument, and tuple of axes as\nkeyword argument ``axis``.\n",
            "metadata": {}
        },
        {
            "id": "d4897098",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy import datasets, ndimage\nfrom scipy.ndimage import vectorized_filter\nimport matplotlib.pyplot as plt\nascent = ndimage.zoom(datasets.ascent(), 0.5).astype(np.float16)\nascent[::16, ::16] = np.nan\nresult = vectorized_filter(ascent, function=np.nanmedian, size=4)",
            "outputs": []
        },
        {
            "id": "9329857a",
            "cell_type": "markdown",
            "source": "Plot the original and filtered images.\n",
            "metadata": {}
        },
        {
            "id": "6756681a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "fig = plt.figure()\nplt.gray()  # show the filtered result in grayscale\nax1 = fig.add_subplot(121)  # left side\nax2 = fig.add_subplot(122)  # right side\nax1.imshow(ascent)\nax2.imshow(result)\nfig.tight_layout()\nplt.show()",
            "outputs": []
        },
        {
            "id": "978c5d11",
            "cell_type": "markdown",
            "source": "Another need satisfied by `vectorized_filter` is to perform multi-output\nfilters. For instance, suppose we wish to filter an image according to the 25th\nand 75th percentiles in addition to the median. We could perform the three\nfilters separately.\n",
            "metadata": {}
        },
        {
            "id": "90e1a923",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "ascent = ndimage.zoom(datasets.ascent(), 0.5)\ndef get_quantile_fun(p):\n    return lambda x, axis: np.quantile(x, p, axis=axis)\nref1 = vectorized_filter(ascent, get_quantile_fun(0.25), size=4)\nref2 = vectorized_filter(ascent, get_quantile_fun(0.50), size=4)\nref3 = vectorized_filter(ascent, get_quantile_fun(0.75), size=4)\nref = np.stack([ref1, ref2, ref3])",
            "outputs": []
        },
        {
            "id": "3d39051f",
            "cell_type": "markdown",
            "source": "However, `vectorized_filter` also supports filters that return multiple outputs\nas long as `output` is unspecified and `batch_memory` is sufficiently high to\nperform the calculation in a single chunk.\n",
            "metadata": {}
        },
        {
            "id": "e734983c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def quartiles(x, axis):\n    return np.quantile(x, [0.25, 0.50, 0.75], axis=axis)\nres = vectorized_filter(ascent, quartiles, size=4, batch_memory=np.inf)\nnp.all(np.isclose(res, ref))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "np.True_"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "f8a56548",
            "cell_type": "markdown",
            "source": "The placement of the additional dimension(s) corresponding with multiple outputs\nis at the discretion of `function`. `quartiles` happens to prepend one dimension\ncorresponding with the three outputs simply because that is the behavior of\n`np.quantile`:\n",
            "metadata": {}
        },
        {
            "id": "c03ad20d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res.shape == (3,) + ascent.shape",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "032f9538",
            "cell_type": "markdown",
            "source": "If we wished for this dimension to be appended:\n",
            "metadata": {}
        },
        {
            "id": "7a716e1b",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def quartiles(x, axis):\n    return np.moveaxis(np.quantile(x, [0.25, 0.50, 0.75], axis=axis), 0, -1)\nres = vectorized_filter(ascent, quartiles, size=4, batch_memory=np.inf)\nres.shape == ascent.shape + (3,)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "29f7ae7f",
            "cell_type": "markdown",
            "source": "Suppose we wish to implment a \"mode\" filter - a filter that selects the most\nfrequently occuring value within the window. A simple (but rather slow)\napproach is to use `generic_filter` with `scipy.stats.mode`.\n",
            "metadata": {}
        },
        {
            "id": "dddf734d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy import stats\nrng = np.random.default_rng()\ninput = rng.integers(255, size=(50, 50)).astype(np.uint8)\ndef simple_mode(input):\n    return stats.mode(input, axis=None).mode\nref = ndimage.generic_filter(input, simple_mode, size=5)",
            "outputs": []
        },
        {
            "id": "9e7a4825",
            "cell_type": "markdown",
            "source": "If speed is important, `vectorized_filter` can take advantage of the performance\nbenefit of a vectorized callable.\n",
            "metadata": {}
        },
        {
            "id": "2f7a6499",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def vectorized_mode(x, axis=(-1,)):\n    n_axes = 1 if np.isscalar(axis) else len(axis)\n    x = np.moveaxis(x, axis, tuple(range(-n_axes, 0)))\n    x = np.reshape(x, x.shape[:-n_axes] + (-1,))\n    y = np.sort(x, axis=-1)\n    i = np.concatenate([np.ones(y.shape[:-1] + (1,), dtype=bool),\n                        y[..., :-1] != y[..., 1:]], axis=-1)\n    indices = np.arange(y.size)[i.ravel()]\n    counts = np.diff(indices, append=y.size)\n    counts = np.reshape(np.repeat(counts, counts), y.shape)\n    k = np.argmax(counts, axis=-1, keepdims=True)\n    return np.take_along_axis(y, k, axis=-1)[..., 0]\nres = vectorized_filter(input, vectorized_mode, size=5)\nnp.all(res == ref)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "np.True_"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "2ae57442",
            "cell_type": "markdown",
            "source": "Depending on the machine, the `vectorized_filter` version may be as much as\n100x faster.\n",
            "metadata": {}
        }
    ]
}