"""
Functions for geometrical image transformation and warping.
"""
import sys
import time
import numbers
from functools import wraps
from collections import namedtuple

import numpy as np
import scipy.ndimage
from scipy.signal import convolve2d

from sunpy import log
from sunpy.util.exceptions import warn_deprecated, warn_user

__all__ = ['add_rotation_function', 'affine_transform']


def affine_transform(image, rmatrix, order=3, scale=1.0, image_center=None,
                     recenter=False, missing=np.nan, use_scipy=None, *, method='scipy', clip=True):
    """
    Rotates, shifts and scales an image.

    This function supports NaN values in the input image and supports using NaN for
    pixels in the output image that are beyond the extent of the input image.
    Handling NaN values in the input image requires additional computation time.

    Parameters
    ----------
    image : `numpy.ndarray`
        2D image to be rotated.
    rmatrix : `numpy.ndarray` that is 2x2
        Linear transformation rotation matrix.
    order : `int` 0-5, optional
        Interpolation order to be used, defaults to 3.  The precise meaning depends
        on the rotation method specifed by ``method``.
    scale : `float`
        A scale factor for the image with the default being no scaling.
    image_center : tuple, optional
        The point in the image to rotate around (axis of rotation).
        Defaults to the center of the array.
    recenter : `bool` or array-like, optional
        Move the axis of rotation to the center of the array or recenter coords.
        Defaults to `False`.
    missing : `float`, optional
        The value to use for pixels in the output image that are beyond the extent
        of the input image.
        Defaults to `numpy.nan`.
    method : {{{rotation_function_names}}}, optional
        Rotation function to use.
        Defaults to ``'scipy'``.
    clip : `bool`, optional
        If `True`, clips the pixel values of the output image to the range of the
        input image (including the value of ``missing``, if used).
        Defaults to `True`.

    Returns
    -------
    `numpy.ndarray`:
        New rotated, scaled and translated image.

    Notes
    -----
    For each NaN pixel in the input image, one or more pixels in the output image
    will be set to NaN, with the size of the pixel region affected depending on the
    interpolation order.  All currently implemented rotation methods require a
    convolution step to handle image NaNs.  This convolution normally uses
    :func:`scipy.signal.convolve2d`, but if `OpenCV <https://opencv.org>`__ is
    installed, the faster |cv2_filter2D|_ is used instead.

    See :func:`~sunpy.image.transform.add_rotation_function` for how to add a
    different rotation function.
    """
    rmatrix = rmatrix / scale
    array_center = (np.array(image.shape)[::-1] - 1) / 2.0

    # Make sure the image center is an array and is where it's supposed to be
    if image_center is not None:
        image_center = np.asanyarray(image_center)
    else:
        image_center = array_center

    # Determine center of rotation based on use (or not) of the recenter keyword
    if recenter:
        rot_center = array_center
    else:
        rot_center = image_center

    displacement = np.dot(rmatrix, rot_center)
    shift = image_center - displacement

    # While `use_scipy` is still supported, we have to check which method to actually use
    method = _get_transform_method(method, use_scipy)

    # Transform the image using the appropriate function
    rotated_image = _rotation_registry[method].function(image, rmatrix, shift, order, missing, clip)

    return rotated_image


def _get_transform_method(method, use_scipy):
    # This is re-used in affine_transform and GenericMap.rotate
    if method not in _rotation_registry:
        raise ValueError(f'Method {method} not in supported methods: '
                         f'{_rotation_registry.keys()}')

    if use_scipy is not None:
        warn_deprecated("The 'use_scipy' argument is deprecated. "
                        "Specify the rotation method to the 'method' "
                        "keyword argument instead.")
        if use_scipy is True and method != 'scipy':
            warn_user(f"Using scipy instead of {method} for rotation.")
            method = 'scipy'

    if method == 'scikit-image':
        try:
            import skimage  # NoQA
        except ImportError:
            raise ImportError("scikit-image must be installed to be usable for rotation.")

    return method


def add_rotation_function(name, *, allowed_orders,
                          handles_clipping, handles_image_nans, handles_nan_missing):
    """
    Decorator to add a rotation function to the registry of selectable
    implementations.

    Each registered rotation function becomes a selectable option for
    :func:`sunpy.image.transform.affine_transform` and
    :meth:`sunpy.map.GenericMap.rotate`.  Those two routines are required to handle
    clipping the output image, NaNs in the input image, and NaN as the value to use
    for pixels in the output image that are beyond the extent of the input image.
    If the supplied rotation function cannot provide one or more of these capabilities,
    the decorator is able to provide them instead.

    The decorator requires the parameters listed under ``Parameters``.  The decorated
    rotation function must accept the parameters listed under ``Other Parameters``
    in that order and return the rotated image.

    Parameters
    ----------
    name : `str`
        The name that will be used to select the rotation function
    allowed_orders : `set`
        The allowed values for the ``order`` parameter.
    handles_clipping : `bool`
        Specifies whether the rotation function will internally perform clipping.
        If ``False``, the rotation function will always receive ``False`` for the
        ``clip`` input parameter.
    handles_image_nans : `bool`
        Specifies whether the rotation function will internally handle NaNs in the
        input image.  If ``False``, the rotation function is guaranteed to be
        provided an image without any NaNs.
    handles_nan_missing : `bool`
        Specifies whether the rotation function will internally handle NaN as the
        ``missing`` value.  If ``False``, the rotation function will never receive
        NaN, but instead receive a value in the input range of the image.

    Other Parameters
    ----------------
    image : `numpy.ndarray`
        The image, which could be integers or floats
    matrix : `numpy.ndarray` that is 2x2
        The linear transformation matrix (e.g., rotation+scale+skew)
    shift : 2-element `numpy.ndarray`
        The translational shift to apply to the image in each axis
    order : `int`
        The numerical parameter that controls the degree of interpolation
    missing : `float`
        The value to use for outside the bounds of the original image
    clip : `bool`
        Whether to clip the output image to the range of the input image

    Notes
    -----
    The docstring of the rotation function should be a bulleted list of notes
    specific to the rotation function.  It will be appended to ``Notes`` section of
    the docstring for :func:`~sunpy.image.transform.affine_transform`.

    The rotation function is supplied the input image directly, so the function
    should not modify the image in place.

    Setting any of the ``handles_*`` parameters to ``False`` means that computation
    will be performed to modify the image returned by the rotation function before
    it is returned to :func:`~sunpy.image.transform.affine_transform`.

    If the decorator is handling image NaNs on behalf of the rotation function
    (i.e., ``handles_image_nans=False``), pixels in the output image will be set to
    NaN if they are within a certain neighborhood size that depends on the ``order``
    parameter.  This step requires an additional image convolution, which might be
    avoidable if the rotation function were able to internally handle image NaNs.
    This convolution normally uses :func:`scipy.signal.convolve2d`, but if
    `OpenCV <https://opencv.org>`__ is installed, the faster |cv2_filter2D|_ is
    used instead.
    """
    def decorator(rotation_function):
        @wraps(rotation_function)
        def wrapper(image, matrix, shift, order, missing, clip):
            if order not in allowed_orders:
                raise ValueError(f"{order} is one of the allowed orders for method '{name}': "
                                 f"{set(allowed_orders)}")

            clip_to_use = clip if handles_clipping else False

            # Check if missing is NaN and needs to be externally handled
            needs_missing_handling = not handles_nan_missing and np.isnan(missing)

            # Check if there are any NaNs in the image that need to be externally handled
            if not handles_image_nans:
                isnan = np.isnan(image)
                needs_nan_handling = np.any(isnan)
            else:
                needs_nan_handling = False

            # If either is needed, change the NaNs to the median of the input image
            if needs_missing_handling or needs_nan_handling:
                substitute = np.nanmedian(image)
            missing_to_use = substitute if needs_missing_handling else missing
            image_to_use = np.nan_to_num(image, nan=substitute) if needs_nan_handling else image

            t = time.perf_counter()
            rotated_image = rotation_function(image_to_use, matrix, shift, order,
                                              missing_to_use, clip_to_use)
            log.debug(f"{name} rotating image: {time.perf_counter() - t:.3f} s")

            # If needed, restore the NaNs
            if needs_nan_handling:
                # Use a convolution to find all pixels that are appreciably affected by NaNs
                # The kernel size depends on the interpolation order, but are empirically defined
                # because a given pixel can affect every other pixel under spline interpolation
                sizes = [1, 1, 5, 5, 7, 7]

                t = time.perf_counter()
                try:
                    # If OpenCV is installed, its convolution function is much faster
                    import cv2
                    expanded_nans = cv2.filter2D(isnan.astype(float), -1,
                                                 np.ones((sizes[order], sizes[order])),
                                                 borderType=cv2.BORDER_CONSTANT)
                except ImportError:
                    expanded_nans = convolve2d(isnan.astype(float),
                                               np.ones((sizes[order], sizes[order])),
                                               mode='same')
                log.debug(f"{name} expanding image NaNs: {time.perf_counter() - t:.3f} s")

                t = time.perf_counter()
                rotated_nans = rotation_function(expanded_nans, matrix, shift, order=min(order, 1),
                                                 missing=0, clip=False)
                rotated_image[rotated_nans > 0] = np.nan
                log.debug(f"{name} rotating image NaNs: {time.perf_counter() - t:.3f} s")

            if needs_missing_handling:
                t = time.perf_counter()
                # Rotate a constant image to determine where to apply the NaNs for `missing`
                constant = rotation_function(np.ones_like(image), matrix, shift, order=1,
                                             missing=0, clip=False)
                rotated_image[constant < 1] = np.nan
                log.debug(f"{name} applying NaN missing: {time.perf_counter() - t:.3f} s")

            if not handles_clipping and clip and not np.all(np.isnan(rotated_image)):
                t = time.perf_counter()
                # Clip the image to the input range
                if np.isnan(missing):
                    # If `missing` is NaN, clipping to the input range is straightforward
                    rotated_image.clip(np.nanmin(image), np.nanmax(image), out=rotated_image)
                else:
                    # Otherwise, check if `missing` should be considered part of the input range
                    lower = np.nanmin([np.max([missing, np.nanmin(rotated_image)]), np.nanmin(image)])
                    upper = np.nanmax([np.min([missing, np.nanmax(rotated_image)]), np.nanmax(image)])
                    rotated_image.clip(lower, upper, out=rotated_image)
                log.debug(f"{name} clipping image: {time.perf_counter() - t:.3f} s")

            return rotated_image

        _rotation_registry[name] = _rotation_method(function=wrapper,
                                                    allowed_orders=set(allowed_orders))

        # Add the docstring of the rotation function to the docstring of affine_transform
        affine_transform.__doc__ += (f"\n    **Specific notes for the '{name}' rotation method:**"
                                     f"\n{wrapper.__doc__}")

        return wrapper
    return decorator


_rotation_method = namedtuple('_rotation_method', ['function', 'allowed_orders'])
_rotation_registry = {}


@add_rotation_function("scipy", allowed_orders=range(6),
                       handles_clipping=False, handles_image_nans=False, handles_nan_missing=True)
def _rotation_scipy(image, matrix, shift, order, missing, clip):
    """
    * Rotates using :func:`scipy.ndimage.affine_transform`
    * The ``order`` parameter is the order of the spline interpolation, and ranges
      from 0 to 5.
    * The ``mode`` parameter for :func:`~scipy.ndimage.affine_transform` is fixed to
      be ``'constant'``
    """
    rotated_image = scipy.ndimage.affine_transform(image.T, matrix, offset=shift, order=order,
                                                   mode='constant', cval=missing).T

    return rotated_image


@add_rotation_function("scikit-image", allowed_orders=range(6),
                       handles_clipping=False, handles_image_nans=False, handles_nan_missing=False)
def _rotation_skimage(image, matrix, shift, order, missing, clip):
    """
    * Rotates using :func:`skimage.transform.warp`
    * The ``order`` parameter selects from the following interpolation algorithms:

      * 0: nearest-neighbor
      * 1: bi-linear
      * 2: bi-quadratic
      * 3: bi-cubic
      * 4: bi-quartic
      * 5: bi-quintic

    * The implementation for higher orders of interpolation means that the pixels
      in the output image that are beyond the extent of the input image may not have
      exactly the value of the ``missing`` parameter.
    * An input image with byte ordering that does not match the native byte order of
      the system (e.g., big-endian values on a little-endian system) will be
      copied and byte-swapped prior to rotation.
    * An input image with integer data is cast to floats prior to passing to
      :func:`~skimage.transform.warp`.  The output image can be re-cast using
      :meth:`numpy.ndarray.astype` if desired.
    * Does not let :func:`~skimage.transform.warp` handle clipping due to
      inconsistent handling across interpolation orders
    * Does not let :func:`~skimage.transform.warp` handle image NaNs because they
      are not handled properly for some interpolation orders
    * Does not pass NaN as ``missing`` to :func:`~skimage.transform.warp` due to
      inconsistent handling across interpolation orders
    """
    import skimage.transform

    # Make the rotation matrix 3x3 to include translation of the image
    skmatrix = np.zeros((3, 3))
    skmatrix[:2, :2] = matrix
    skmatrix[2, 2] = 1.0
    skmatrix[:2, 2] = shift
    tform = skimage.transform.AffineTransform(skmatrix)

    if issubclass(image.dtype.type, numbers.Integral):
        warn_user("Integer input data has been cast to float64.")
        adjusted_image = image.astype(np.float64)
    else:
        adjusted_image = image.copy()

    # Scale image to range [0, 1]
    im_min = np.nanmin([missing, np.min(adjusted_image)])
    adjusted_image -= im_min
    adjusted_missing = missing - im_min
    im_max = np.nanmax([adjusted_missing, np.max(adjusted_image)])
    if im_max > 0:
        adjusted_image /= im_max
        adjusted_missing /= im_max

    # Swap the byte order if it is non-native (e.g., big-endian on a little-endian system)
    if adjusted_image.dtype.byteorder == ('>' if sys.byteorder == 'little' else '<'):
        adjusted_image = adjusted_image.byteswap().newbyteorder()

    # Be aware that even though mode is set to 'constant', when skimage 0.19 calls scipy,
    # it specifies the scipy mode to be 'grid-constant' rather than 'constant'
    rotated_image = skimage.transform.warp(adjusted_image, tform, order=order,
                                           mode='constant', cval=adjusted_missing, clip=clip)

    # Convert the image back to its original range
    if im_max > 0:
        rotated_image *= im_max
    rotated_image += im_min

    return rotated_image


@add_rotation_function("opencv", allowed_orders={0, 1, 3},
                       handles_clipping=False, handles_image_nans=False, handles_nan_missing=False)
def _rotation_cv2(image, matrix, shift, order, missing, clip):
    """
    * Rotates using |cv2_warpAffine|_ from `OpenCV <https://opencv.org>`__
    * The ``order`` parameter selects from the following interpolation algorithms:

      * 0: nearest-neighbor interpolation
      * 1: bilinear interpolation
      * 3: bicubic interpolation

    * An input image with byte ordering that does not match the native byte order of
      the system (e.g., big-endian values on a little-endian system) will be
      copied and byte-swapped prior to rotation.
    * An input image with integer data is cast to floats prior to passing to
      |cv2_warpAffine|_.  The output image can be re-cast using
      :meth:`numpy.ndarray.astype` if desired.
    """
    try:
        import cv2
    except ImportError:
        raise ImportError("The opencv-python package is required to use this rotation method.")

    _CV_ORDER_FLAGS = {
        0: cv2.INTER_NEAREST,
        1: cv2.INTER_LINEAR,
        3: cv2.INTER_CUBIC,
    }
    order_to_use = _CV_ORDER_FLAGS[order]

    if issubclass(image.dtype.type, numbers.Integral):
        warn_user("Integer input data has been cast to float64.")
        adjusted_image = image.astype(np.float64)
    else:
        adjusted_image = image.copy()

    trans = np.concatenate([matrix, shift[:, np.newaxis]], axis=1)

    # Swap the byte order if it is non-native (e.g., big-endian on a little-endian system)
    if adjusted_image.dtype.byteorder == ('>' if sys.byteorder == 'little' else '<'):
        adjusted_image = adjusted_image.byteswap().newbyteorder()

    # missing must be a Python float, not a NumPy float
    missing = float(missing)

    return cv2.warpAffine(adjusted_image, trans, np.flip(adjusted_image.shape),
                          flags=order_to_use | cv2.WARP_INVERSE_MAP,
                          borderMode=cv2.BORDER_CONSTANT, borderValue=missing)


# Generate the string with allowable rotation-function names for use in docstrings
_rotation_function_names = ", ".join([f"``'{name}'``" for name in _rotation_registry])
# Insert into the docstring for affine_transform.  We cannot use the add_common_docstring decorator
# due to what would be a circular loop in definitions.
affine_transform.__doc__ = affine_transform.__doc__.format(rotation_function_names=_rotation_function_names)

"""
Functions for geometrical image transformation and warping.
"""
import sys
import time
import numbers
from functools import wraps
from collections import namedtuple

import numpy as np
import scipy.ndimage
from scipy.signal import convolve2d

from sunpy import log
from sunpy.util.exceptions import warn_deprecated, warn_user

__all__ = ['add_rotation_function', 'affine_transform']


def affine_transform(image, rmatrix, order=3, scale=1.0, image_center=None,
                     recenter=False, missing=np.nan, use_scipy=None, *, method='scipy', clip=True):
    """
    Rotates, shifts and scales an image.

    This function supports NaN values in the input image and supports using NaN for
    pixels in the output image that are beyond the extent of the input image.
    Handling NaN values in the input image requires additional computation time.

    Parameters
    ----------
    image : `numpy.ndarray`
        2D image to be rotated.
    rmatrix : `numpy.ndarray` that is 2x2
        Linear transformation rotation matrix.
    order : `int` 0-5, optional
        Interpolation order to be used, defaults to 3.  The precise meaning depends
        on the rotation method specifed by ``method``.
    scale : `float`
        A scale factor for the image with the default being no scaling.
    image_center : tuple, optional
        The point in the image to rotate around (axis of rotation).
        Defaults to the center of the array.
    recenter : `bool` or array-like, optional
        Move the axis of rotation to the center of the array or recenter coords.
        Defaults to `False`.
    missing : `float`, optional
        The value to use for pixels in the output image that are beyond the extent
        of the input image.
        Defaults to `numpy.nan`.
    method : {{{rotation_function_names}}}, optional
        Rotation function to use.
        Defaults to ``'scipy'``.
    clip : `bool`, optional
        If `True`, clips the pixel values of the output image to the range of the
        input image (including the value of ``missing``, if used).
        Defaults to `True`.

    Returns
    -------
    `numpy.ndarray`:
        New rotated, scaled and translated image.

    Notes
    -----
    For each NaN pixel in the input image, one or more pixels in the output image
    will be set to NaN, with the size of the pixel region affected depending on the
    interpolation order.  All currently implemented rotation methods require a
    convolution step to handle image NaNs.  This convolution normally uses
    :func:`scipy.signal.convolve2d`, but if `OpenCV <https://opencv.org>`__ is
    installed, the faster |cv2_filter2D|_ is used instead.

    See :func:`~sunpy.image.transform.add_rotation_function` for how to add a
    different rotation function.
    """
    rmatrix = rmatrix / scale
    array_center = (np.array(image.shape)[::-1] - 1) / 2.0

    # Make sure the image center is an array and is where it's supposed to be
    if image_center is not None:
        image_center = np.asanyarray(image_center)
    else:
        image_center = array_center

    # Determine center of rotation based on use (or not) of the recenter keyword
    if recenter:
        rot_center = array_center
    else:
        rot_center = image_center

    displacement = np.dot(rmatrix, rot_center)
    shift = image_center - displacement

    # While `use_scipy` is still supported, we have to check which method to actually use
    method = _get_transform_method(method, use_scipy)

    # Transform the image using the appropriate function
    rotated_image = _rotation_registry[method].function(image, rmatrix, shift, order, missing, clip)

    return rotated_image


def _get_transform_method(method, use_scipy):
    # This is re-used in affine_transform and GenericMap.rotate
    if method not in _rotation_registry:
        raise ValueError(f'Method {method} not in supported methods: '
                         f'{_rotation_registry.keys()}')

    if use_scipy is not None:
        warn_deprecated("The 'use_scipy' argument is deprecated. "
                        "Specify the rotation method to the 'method' "
                        "keyword argument instead.")
        if use_scipy is True and method != 'scipy':
            warn_user(f"Using scipy instead of {method} for rotation.")
            method = 'scipy'

    if method == 'scikit-image':
        try:
            import skimage  # NoQA
        except ImportError:
            raise ImportError("scikit-image must be installed to be usable for rotation.")

    return method


def add_rotation_function(name, *, allowed_orders,
                          handles_clipping, handles_image_nans, handles_nan_missing):
    """
    Decorator to add a rotation function to the registry of selectable
    implementations.

    Each registered rotation function becomes a selectable option for
    :func:`sunpy.image.transform.affine_transform` and
    :meth:`sunpy.map.GenericMap.rotate`.  Those two routines are required to handle
    clipping the output image, NaNs in the input image, and NaN as the value to use
    for pixels in the output image that are beyond the extent of the input image.
    If the supplied rotation function cannot provide one or more of these capabilities,
    the decorator is able to provide them instead.

    The decorator requires the parameters listed under ``Parameters``.  The decorated
    rotation function must accept the parameters listed under ``Other Parameters``
    in that order and return the rotated image.

    Parameters
    ----------
    name : `str`
        The name that will be used to select the rotation function
    allowed_orders : `set`
        The allowed values for the ``order`` parameter.
    handles_clipping : `bool`
        Specifies whether the rotation function will internally perform clipping.
        If ``False``, the rotation function will always receive ``False`` for the
        ``clip`` input parameter.
    handles_image_nans : `bool`
        Specifies whether the rotation function will internally handle NaNs in the
        input image.  If ``False``, the rotation function is guaranteed to be
        provided an image without any NaNs.
    handles_nan_missing : `bool`
        Specifies whether the rotation function will internally handle NaN as the
        ``missing`` value.  If ``False``, the rotation function will never receive
        NaN, but instead receive a value in the input range of the image.

    Other Parameters
    ----------------
    image : `numpy.ndarray`
        The image, which could be integers or floats
    matrix : `numpy.ndarray` that is 2x2
        The linear transformation matrix (e.g., rotation+scale+skew)
    shift : 2-element `numpy.ndarray`
        The translational shift to apply to the image in each axis
    order : `int`
        The numerical parameter that controls the degree of interpolation
    missing : `float`
        The value to use for outside the bounds of the original image
    clip : `bool`
        Whether to clip the output image to the range of the input image

    Notes
    -----
    The docstring of the rotation function should be a bulleted list of notes
    specific to the rotation function.  It will be appended to ``Notes`` section of
    the docstring for :func:`~sunpy.image.transform.affine_transform`.

    The rotation function is supplied the input image directly, so the function
    should not modify the image in place.

    Setting any of the ``handles_*`` parameters to ``False`` means that computation
    will be performed to modify the image returned by the rotation function before
    it is returned to :func:`~sunpy.image.transform.affine_transform`.

    If the decorator is handling image NaNs on behalf of the rotation function
    (i.e., ``handles_image_nans=False``), pixels in the output image will be set to
    NaN if they are within a certain neighborhood size that depends on the ``order``
    parameter.  This step requires an additional image convolution, which might be
    avoidable if the rotation function were able to internally handle image NaNs.
    This convolution normally uses :func:`scipy.signal.convolve2d`, but if
    `OpenCV <https://opencv.org>`__ is installed, the faster |cv2_filter2D|_ is
    used instead.
    """
    def decorator(rotation_function):
        @wraps(rotation_function)
        def wrapper(image, matrix, shift, order, missing, clip):
            if order not in allowed_orders:
                raise ValueError(f"{order} is one of the allowed orders for method '{name}': "
                                 f"{set(allowed_orders)}")

            clip_to_use = clip if handles_clipping else False

            # Check if missing is NaN and needs to be externally handled
            needs_missing_handling = not handles_nan_missing and np.isnan(missing)

            # Check if there are any NaNs in the image that need to be externally handled
            if not handles_image_nans:
                isnan = np.isnan(image)
                needs_nan_handling = np.any(isnan)
            else:
                needs_nan_handling = False

            # If either is needed, change the NaNs to the median of the input image
            if needs_missing_handling or needs_nan_handling:
                substitute = np.nanmedian(image)
            missing_to_use = substitute if needs_missing_handling else missing
            image_to_use = np.nan_to_num(image, nan=substitute) if needs_nan_handling else image

            t = time.perf_counter()
            rotated_image = rotation_function(image_to_use, matrix, shift, order,
                                              missing_to_use, clip_to_use)
            log.debug(f"{name} rotating image: {time.perf_counter() - t:.3f} s")

            # If needed, restore the NaNs
            if needs_nan_handling:
                # Use a convolution to find all pixels that are appreciably affected by NaNs
                # The kernel size depends on the interpolation order, but are empirically defined
                # because a given pixel can affect every other pixel under spline interpolation
                sizes = [1, 1, 5, 5, 7, 7]

                t = time.perf_counter()
                try:
                    # If OpenCV is installed, its convolution function is much faster
                    import cv2
                    expanded_nans = cv2.filter2D(isnan.astype(float), -1,
                                                 np.ones((sizes[order], sizes[order])),
                                                 borderType=cv2.BORDER_CONSTANT)
                except ImportError:
                    expanded_nans = convolve2d(isnan.astype(float),
                                               np.ones((sizes[order], sizes[order])),
                                               mode='same')
                log.debug(f"{name} expanding image NaNs: {time.perf_counter() - t:.3f} s")

                t = time.perf_counter()
                rotated_nans = rotation_function(expanded_nans, matrix, shift, order=min(order, 1),
                                                 missing=0, clip=False)
                rotated_image[rotated_nans > 0] = np.nan
                log.debug(f"{name} rotating image NaNs: {time.perf_counter() - t:.3f} s")

            if needs_missing_handling:
                t = time.perf_counter()
                # Rotate a constant image to determine where to apply the NaNs for `missing`
                constant = rotation_function(np.ones_like(image), matrix, shift, order=1,
                                             missing=0, clip=False)
                rotated_image[constant < 1] = np.nan
                log.debug(f"{name} applying NaN missing: {time.perf_counter() - t:.3f} s")

            if not handles_clipping and clip and not np.all(np.isnan(rotated_image)):
                t = time.perf_counter()
                # Clip the image to the input range
                if np.isnan(missing):
                    # If `missing` is NaN, clipping to the input range is straightforward
                    rotated_image.clip(np.nanmin(image), np.nanmax(image), out=rotated_image)
                else:
                    # Otherwise, check if `missing` should be considered part of the input range
                    lower = np.nanmin([np.max([missing, np.nanmin(rotated_image)]), np.nanmin(image)])
                    upper = np.nanmax([np.min([missing, np.nanmax(rotated_image)]), np.nanmax(image)])
                    rotated_image.clip(lower, upper, out=rotated_image)
                log.debug(f"{name} clipping image: {time.perf_counter() - t:.3f} s")

            return rotated_image

        _rotation_registry[name] = _rotation_method(function=wrapper,
                                                    allowed_orders=set(allowed_orders))

        # Add the docstring of the rotation function to the docstring of affine_transform
        affine_transform.__doc__ += (f"\n    **Specific notes for the '{name}' rotation method:**"
                                     f"\n{wrapper.__doc__}")

        return wrapper
    return decorator


_rotation_method = namedtuple('_rotation_method', ['function', 'allowed_orders'])
_rotation_registry = {}


@add_rotation_function("scipy", allowed_orders=range(6),
                       handles_clipping=False, handles_image_nans=False, handles_nan_missing=True)
def _rotation_scipy(image, matrix, shift, order, missing, clip):
    """
    * Rotates using :func:`scipy.ndimage.affine_transform`
    * The ``order`` parameter is the order of the spline interpolation, and ranges
      from 0 to 5.
    * The ``mode`` parameter for :func:`~scipy.ndimage.affine_transform` is fixed to
      be ``'constant'``
    """
    rotated_image = scipy.ndimage.affine_transform(image.T, matrix, offset=shift, order=order,
                                                   mode='constant', cval=missing).T

    return rotated_image


@add_rotation_function("scikit-image", allowed_orders=range(6),
                       handles_clipping=False, handles_image_nans=False, handles_nan_missing=False)
def _rotation_skimage(image, matrix, shift, order, missing, clip):
    """
    * Rotates using :func:`skimage.transform.warp`
    * The ``order`` parameter selects from the following interpolation algorithms:

      * 0: nearest-neighbor
      * 1: bi-linear
      * 2: bi-quadratic
      * 3: bi-cubic
      * 4: bi-quartic
      * 5: bi-quintic

    * The implementation for higher orders of interpolation means that the pixels
      in the output image that are beyond the extent of the input image may not have
      exactly the value of the ``missing`` parameter.
    * An input image with byte ordering that does not match the native byte order of
      the system (e.g., big-endian values on a little-endian system) will be
      copied and byte-swapped prior to rotation.
    * An input image with integer data is cast to floats prior to passing to
      :func:`~skimage.transform.warp`.  The output image can be re-cast using
      :meth:`numpy.ndarray.astype` if desired.
    * Does not let :func:`~skimage.transform.warp` handle clipping due to
      inconsistent handling across interpolation orders
    * Does not let :func:`~skimage.transform.warp` handle image NaNs because they
      are not handled properly for some interpolation orders
    * Does not pass NaN as ``missing`` to :func:`~skimage.transform.warp` due to
      inconsistent handling across interpolation orders
    """
    import skimage.transform

    # Make the rotation matrix 3x3 to include translation of the image
    skmatrix = np.zeros((3, 3))
    skmatrix[:2, :2] = matrix
    skmatrix[2, 2] = 1.0
    skmatrix[:2, 2] = shift
    tform = skimage.transform.AffineTransform(skmatrix)

    if issubclass(image.dtype.type, numbers.Integral):
        warn_user("Integer input data has been cast to float64.")
        adjusted_image = image.astype(np.float64)
    else:
        adjusted_image = image.copy()

    # Scale image to range [0, 1]
    im_min = np.nanmin([missing, np.min(adjusted_image)])
    adjusted_image -= im_min
    adjusted_missing = missing - im_min
    im_max = np.nanmax([adjusted_missing, np.max(adjusted_image)])
    if im_max > 0:
        adjusted_image /= im_max
        adjusted_missing /= im_max

    # Swap the byte order if it is non-native (e.g., big-endian on a little-endian system)
    if adjusted_image.dtype.byteorder == ('>' if sys.byteorder == 'little' else '<'):
        adjusted_image = adjusted_image.byteswap().newbyteorder()

    # Be aware that even though mode is set to 'constant', when skimage 0.19 calls scipy,
    # it specifies the scipy mode to be 'grid-constant' rather than 'constant'
    rotated_image = skimage.transform.warp(adjusted_image, tform, order=order,
                                           mode='constant', cval=adjusted_missing, clip=clip)

    # Convert the image back to its original range
    if im_max > 0:
        rotated_image *= im_max
    rotated_image += im_min

    return rotated_image


@add_rotation_function("opencv", allowed_orders={0, 1, 3},
                       handles_clipping=False, handles_image_nans=False, handles_nan_missing=False)
def _rotation_cv2(image, matrix, shift, order, missing, clip):
    """
    * Rotates using |cv2_warpAffine|_ from `OpenCV <https://opencv.org>`__
    * The ``order`` parameter selects from the following interpolation algorithms:

      * 0: nearest-neighbor interpolation
      * 1: bilinear interpolation
      * 3: bicubic interpolation

    * An input image with byte ordering that does not match the native byte order of
      the system (e.g., big-endian values on a little-endian system) will be
      copied and byte-swapped prior to rotation.
    * An input image with integer data is cast to floats prior to passing to
      |cv2_warpAffine|_.  The output image can be re-cast using
      :meth:`numpy.ndarray.astype` if desired.
    """
    try:
        import cv2
    except ImportError:
        raise ImportError("The opencv-python package is required to use this rotation method.")

    _CV_ORDER_FLAGS = {
        0: cv2.INTER_NEAREST,
        1: cv2.INTER_LINEAR,
        3: cv2.INTER_CUBIC,
    }
    order_to_use = _CV_ORDER_FLAGS[order]

    if issubclass(image.dtype.type, numbers.Integral):
        warn_user("Integer input data has been cast to float64.")
        adjusted_image = image.astype(np.float64)
    else:
        adjusted_image = image.copy()

    trans = np.concatenate([matrix, shift[:, np.newaxis]], axis=1)

    # Swap the byte order if it is non-native (e.g., big-endian on a little-endian system)
    if adjusted_image.dtype.byteorder == ('>' if sys.byteorder == 'little' else '<'):
        adjusted_image = adjusted_image.byteswap().newbyteorder()

    # missing must be a Python float, not a NumPy float
    missing = float(missing)

    return cv2.warpAffine(adjusted_image, trans, np.flip(adjusted_image.shape),
                          flags=order_to_use | cv2.WARP_INVERSE_MAP,
                          borderMode=cv2.BORDER_CONSTANT, borderValue=missing)


# Generate the string with allowable rotation-function names for use in docstrings
_rotation_function_names = ", ".join([f"``'{name}'``" for name in _rotation_registry])
# Insert into the docstring for affine_transform.  We cannot use the add_common_docstring decorator
# due to what would be a circular loop in definitions.
affine_transform.__doc__ = affine_transform.__doc__.format(rotation_function_names=_rotation_function_names)

Functions for geometrical image transformation and warping. Rotates, shifts and scales an image.      This function supports NaN values in the input image and supports using NaN for     pixels in the output image that are beyond the extent of the input image.     Handling NaN values in the input image requires additional computation time.      Parameters     ----------     image : `numpy.ndarray`         2D image to be rotated.     rmatrix : `numpy.ndarray` that is 2x2         Linear transformation rotation matrix.     order : `int` 0-5, optional         Interpolation order to be used, defaults to 3.  The precise meaning depends         on the rotation method specifed by ``method``.     scale : `float`         A scale factor for the image with the default being no scaling.     image_center : tuple, optional         The point in the image to rotate around (axis of rotation).         Defaults to the center of the array.     recenter : `bool` or array-like, optional         Move the axis of rotation to the center of the array or recenter coords.         Defaults to `False`.     missing : `float`, optional         The value to use for pixels in the output image that are beyond the extent         of the input image.         Defaults to `numpy.nan`.     method : {{{rotation_function_names}}}, optional         Rotation function to use.         Defaults to ``'scipy'``.     clip : `bool`, optional         If `True`, clips the pixel values of the output image to the range of the         input image (including the value of ``missing``, if used).         Defaults to `True`.      Returns     -------     `numpy.ndarray`:         New rotated, scaled and translated image.      Notes     -----     For each NaN pixel in the input image, one or more pixels in the output image     will be set to NaN, with the size of the pixel region affected depending on the     interpolation order.  All currently implemented rotation methods require a     convolution step to handle image NaNs.  This convolution normally uses     :func:`scipy.signal.convolve2d`, but if `OpenCV <https://opencv.org>`__ is     installed, the faster |cv2_filter2D|_ is used instead.      See :func:`~sunpy.image.transform.add_rotation_function` for how to add a     different rotation function. # Make sure the image center is an array and is where it's supposed to be # Determine center of rotation based on use (or not) of the recenter keyword # While `use_scipy` is still supported, we have to check which method to actually use # Transform the image using the appropriate function # This is re-used in affine_transform and GenericMap.rotate # NoQA Decorator to add a rotation function to the registry of selectable     implementations.      Each registered rotation function becomes a selectable option for     :func:`sunpy.image.transform.affine_transform` and     :meth:`sunpy.map.GenericMap.rotate`.  Those two routines are required to handle     clipping the output image, NaNs in the input image, and NaN as the value to use     for pixels in the output image that are beyond the extent of the input image.     If the supplied rotation function cannot provide one or more of these capabilities,     the decorator is able to provide them instead.      The decorator requires the parameters listed under ``Parameters``.  The decorated     rotation function must accept the parameters listed under ``Other Parameters``     in that order and return the rotated image.      Parameters     ----------     name : `str`         The name that will be used to select the rotation function     allowed_orders : `set`         The allowed values for the ``order`` parameter.     handles_clipping : `bool`         Specifies whether the rotation function will internally perform clipping.         If ``False``, the rotation function will always receive ``False`` for the         ``clip`` input parameter.     handles_image_nans : `bool`         Specifies whether the rotation function will internally handle NaNs in the         input image.  If ``False``, the rotation function is guaranteed to be         provided an image without any NaNs.     handles_nan_missing : `bool`         Specifies whether the rotation function will internally handle NaN as the         ``missing`` value.  If ``False``, the rotation function will never receive         NaN, but instead receive a value in the input range of the image.      Other Parameters     ----------------     image : `numpy.ndarray`         The image, which could be integers or floats     matrix : `numpy.ndarray` that is 2x2         The linear transformation matrix (e.g., rotation+scale+skew)     shift : 2-element `numpy.ndarray`         The translational shift to apply to the image in each axis     order : `int`         The numerical parameter that controls the degree of interpolation     missing : `float`         The value to use for outside the bounds of the original image     clip : `bool`         Whether to clip the output image to the range of the input image      Notes     -----     The docstring of the rotation function should be a bulleted list of notes     specific to the rotation function.  It will be appended to ``Notes`` section of     the docstring for :func:`~sunpy.image.transform.affine_transform`.      The rotation function is supplied the input image directly, so the function     should not modify the image in place.      Setting any of the ``handles_*`` parameters to ``False`` means that computation     will be performed to modify the image returned by the rotation function before     it is returned to :func:`~sunpy.image.transform.affine_transform`.      If the decorator is handling image NaNs on behalf of the rotation function     (i.e., ``handles_image_nans=False``), pixels in the output image will be set to     NaN if they are within a certain neighborhood size that depends on the ``order``     parameter.  This step requires an additional image convolution, which might be     avoidable if the rotation function were able to internally handle image NaNs.     This convolution normally uses :func:`scipy.signal.convolve2d`, but if     `OpenCV <https://opencv.org>`__ is installed, the faster |cv2_filter2D|_ is     used instead. # Check if missing is NaN and needs to be externally handled # Check if there are any NaNs in the image that need to be externally handled # If either is needed, change the NaNs to the median of the input image # If needed, restore the NaNs # Use a convolution to find all pixels that are appreciably affected by NaNs # The kernel size depends on the interpolation order, but are empirically defined # because a given pixel can affect every other pixel under spline interpolation # If OpenCV is installed, its convolution function is much faster # Rotate a constant image to determine where to apply the NaNs for `missing` # Clip the image to the input range # If `missing` is NaN, clipping to the input range is straightforward # Otherwise, check if `missing` should be considered part of the input range # Add the docstring of the rotation function to the docstring of affine_transform * Rotates using :func:`scipy.ndimage.affine_transform`     * The ``order`` parameter is the order of the spline interpolation, and ranges       from 0 to 5.     * The ``mode`` parameter for :func:`~scipy.ndimage.affine_transform` is fixed to       be ``'constant'`` * Rotates using :func:`skimage.transform.warp`     * The ``order`` parameter selects from the following interpolation algorithms:        * 0: nearest-neighbor       * 1: bi-linear       * 2: bi-quadratic       * 3: bi-cubic       * 4: bi-quartic       * 5: bi-quintic      * The implementation for higher orders of interpolation means that the pixels       in the output image that are beyond the extent of the input image may not have       exactly the value of the ``missing`` parameter.     * An input image with byte ordering that does not match the native byte order of       the system (e.g., big-endian values on a little-endian system) will be       copied and byte-swapped prior to rotation.     * An input image with integer data is cast to floats prior to passing to       :func:`~skimage.transform.warp`.  The output image can be re-cast using       :meth:`numpy.ndarray.astype` if desired.     * Does not let :func:`~skimage.transform.warp` handle clipping due to       inconsistent handling across interpolation orders     * Does not let :func:`~skimage.transform.warp` handle image NaNs because they       are not handled properly for some interpolation orders     * Does not pass NaN as ``missing`` to :func:`~skimage.transform.warp` due to       inconsistent handling across interpolation orders # Make the rotation matrix 3x3 to include translation of the image # Scale image to range [0, 1] # Swap the byte order if it is non-native (e.g., big-endian on a little-endian system) # Be aware that even though mode is set to 'constant', when skimage 0.19 calls scipy, # it specifies the scipy mode to be 'grid-constant' rather than 'constant' # Convert the image back to its original range * Rotates using |cv2_warpAffine|_ from `OpenCV <https://opencv.org>`__     * The ``order`` parameter selects from the following interpolation algorithms:        * 0: nearest-neighbor interpolation       * 1: bilinear interpolation       * 3: bicubic interpolation      * An input image with byte ordering that does not match the native byte order of       the system (e.g., big-endian values on a little-endian system) will be       copied and byte-swapped prior to rotation.     * An input image with integer data is cast to floats prior to passing to       |cv2_warpAffine|_.  The output image can be re-cast using       :meth:`numpy.ndarray.astype` if desired. # Swap the byte order if it is non-native (e.g., big-endian on a little-endian system) # missing must be a Python float, not a NumPy float # Generate the string with allowable rotation-function names for use in docstrings # Insert into the docstring for affine_transform.  We cannot use the add_common_docstring decorator # due to what would be a circular loop in definitions.