#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""Various useful functions for creating filters and textures
(e.g. for PatchStim)
"""

# Part of the PsychoPy library
# Copyright (C) 2002-2018 <NAME> (C) 2019-2021 Open Science Tools Ltd.
# Distributed under the terms of the GNU General Public License (GPL).

from __future__ import absolute_import, division, print_function

from past.utils import old_div
import numpy
from numpy.fft import fft2, ifft2, fftshift, ifftshift
from psychopy import logging
try:
    from PIL import Image
except ImportError:
    import Image


def makeGrating(res,
                ori=0.0,  # in degrees
                cycles=1.0,
                phase=0.0,  # in degrees
                gratType="sin",
                contr=1.0):
    """Make an array containing a luminance grating of the specified params

    :Parameters:
        res: integer
            the size of the resulting matrix on both dimensions (e.g 256)
        ori: float or int (default=0.0)
            the orientation of the grating in degrees
        cycles:float or int (default=1.0)
            the number of grating cycles within the array
        phase: float or int (default=0.0)
            the phase of the grating in degrees (NB this differs to most
            PsychoPy phase arguments which use units of fraction of a cycle)
        gratType: 'sin', 'sqr', 'ramp' or 'sinXsin' (default="sin")
            the type of grating to be 'drawn'
        contr: float (default=1.0)
            contrast of the grating

    :Returns:
        a square numpy array of size resXres

    """
    # to prevent the sinusoid ever being exactly at zero (for sqr wave):
    tiny = 0.0000000000001
    ori *= (old_div(-numpy.pi, 180))
    phase *= (old_div(numpy.pi, 180))
    cyclesTwoPi = cycles * 2.0 * numpy.pi
    xrange, yrange = numpy.mgrid[0.0: cyclesTwoPi: old_div(cyclesTwoPi, res),
                                 0.0: cyclesTwoPi: old_div(cyclesTwoPi, res)]

    sin, cos = numpy.sin, numpy.cos
    if gratType == "none":
        res = 2
        intensity = numpy.ones((res, res), float)
    elif gratType == "sin":
        intensity = contr * sin(xrange * sin(ori) + yrange * cos(ori) + phase)
    elif gratType == "ramp":
        intensity = contr * (xrange * cos(ori) +
                             yrange * sin(ori)) / (2 * numpy.pi)
    elif gratType == "sqr":  # square wave (symmetric duty cycle)
        intensity = numpy.where(sin(xrange * sin(ori) + yrange * cos(ori) +
                                    phase + tiny) >= 0, 1, -1)
    elif gratType == "sinXsin":
        intensity = sin(xrange) * sin(yrange)
    else:
        # might be a filename of an image
        try:
            im = Image.open(gratType)
        except Exception:
            logging.error("couldn't find tex...", gratType)
            return
        # todo: opened it, now what?
    return intensity


def maskMatrix(matrix, shape='circle', radius=1.0, center=(0.0, 0.0)):
    """Make and apply a mask to an input matrix (e.g. a grating)

    :Parameters:
         matrix:  a square numpy array
             array to which the mask should be applied
         shape:  'circle','gauss','ramp' (linear gradient from center)
             shape of the mask
         radius:  float
             scale factor to be applied to the mask (circle with radius of
             [1,1] will extend just to the edge of the matrix). Radius can
             be asymmetric, e.g. [1.0,2.0] will be wider than it is tall.
         center:  2x1 tuple or list (default=[0.0,0.0])
             the centre of the mask in the matrix ([1,1] is top-right
             corner, [-1,-1] is bottom-left)
    """
    # NB makeMask now returns a value -1:1
    alphaMask = makeMask(matrix.shape[0], shape, radius,
                         center=(0.0, 0.0), range=[0, 1])
    return matrix * alphaMask


def makeMask(matrixSize, shape='circle', radius=1.0, center=(0.0, 0.0),
             range=(-1, 1), fringeWidth=0.2):
    """Returns a matrix to be used as an alpha mask (circle,gauss,ramp).

    :Parameters:
            matrixSize: integer
                the size of the resulting matrix on both dimensions (e.g 256)
            shape:  'circle','gauss','ramp' (linear gradient from center),
                'raisedCosine' (the edges are blurred by a raised cosine)
                shape of the mask
            radius:  float
                scale factor to be applied to the mask (circle with radius of
                [1,1] will extend just to the edge of the matrix). Radius can
                asymmetric, e.g. [1.0,2.0] will be wider than it is tall.
            center:  2x1 tuple or list (default=[0.0,0.0])
                the centre of the mask in the matrix ([1,1] is top-right
                corner, [-1,-1] is bottom-left)
            fringeWidth: float (0-1)
                The proportion of the raisedCosine that is being blurred.
            range: 2x1 tuple or list (default=[-1,1])
                The minimum and maximum value in the mask matrix
    """
    rad = makeRadialMatrix(matrixSize, center, radius)
    if shape == 'ramp':
        outArray = 1 - rad
    elif shape == 'circle':
        # outArray=numpy.ones(matrixSize,'f')
        outArray = numpy.where(numpy.greater(rad, 1.0), 0.0, 1.0)
    elif shape == 'gauss':
        outArray = makeGauss(rad, mean=0.0, sd=0.33333)
    elif shape == 'raisedCosine':
        hammingLen = 1000  # This affects the 'granularity' of the raised cos
        fringeProportion = fringeWidth  # This one affects the proportion of
        # the stimulus diameter that is devoted to the raised cosine.

        rad = makeRadialMatrix(matrixSize, center, radius)
        outArray = numpy.zeros_like(rad)
        outArray[numpy.where(rad < 1)] = 1
        raisedCosIdx = numpy.where(
            [numpy.logical_and(rad <= 1, rad >= 1 - fringeProportion)])[1:]

        # Make a raised_cos (half a hamming window):
        raisedCos = numpy.hamming(hammingLen)[:hammingLen//2]
        raisedCos -= numpy.min(raisedCos)
        raisedCos /= numpy.max(raisedCos)

        # Measure the distance from the edge - this is your index into the
        # hamming window:
        dFromEdge = numpy.abs((1 - fringeProportion) - rad[raisedCosIdx])
        dFromEdge /= numpy.max(dFromEdge)
        dFromEdge *= numpy.round(hammingLen/2)

        # This is the indices into the hamming (larger for small distances
        # from the edge!):
        portion_idx = (-1 * dFromEdge).astype(int)

        # Apply the raised cos to this portion:
        outArray[raisedCosIdx] = raisedCos[portion_idx]

        # Sometimes there are some remaining artifacts from this process, get
        # rid of them:
        artifact_idx = numpy.where(
            numpy.logical_and(outArray == 0, rad < 0.99))
        outArray[artifact_idx] = 1
        artifact_idx = numpy.where(
            numpy.logical_and(outArray == 1, rad > 0.99))
        outArray[artifact_idx] = 0

    else:
        raise ValueError('Unknown value for shape argument %s' % shape)
    mag = range[1] - range[0]
    offset = range[0]
    return outArray * mag + offset


def makeRadialMatrix(matrixSize, center=(0.0, 0.0), radius=1.0):
    """Generate a square matrix where each element val is
    its distance from the centre of the matrix

    :Parameters:
        matrixSize: integer
            the size of the resulting matrix on both dimensions (e.g 256)
        radius:  float
            scale factor to be applied to the mask (circle with radius of
            [1,1] will extend just to the edge of the matrix). Radius can
            be asymmetric, e.g. [1.0,2.0] will be wider than it is tall.
        center:  2x1 tuple or list (default=[0.0,0.0])
            the centre of the mask in the matrix ([1,1] is top-right
            corner, [-1,-1] is bottom-left)
    """
    if type(radius) in [int, float]:
        radius = [radius, radius]

    # NB need to add one step length because
    yy, xx = numpy.mgrid[0:matrixSize, 0:matrixSize]
    xx = ((1.0 - 2.0 / matrixSize * xx) + center[0]) / radius[0]
    yy = ((1.0 - 2.0 / matrixSize * yy) + center[1]) / radius[1]
    rad = numpy.sqrt(numpy.power(xx, 2) + numpy.power(yy, 2))
    return rad


def makeGauss(x, mean=0.0, sd=1.0, gain=1.0, base=0.0):
    """
    Return the gaussian distribution for a given set of x-vals

   :Parameters:
        mean: float
            the centre of the distribution
        sd: float
            the width of the distribution
        gain: float
            the height of the distribution
        base: float
            an offset added to the result

    """
    simpleGauss = numpy.exp((-numpy.power(mean - x, 2) / (2 * sd**2)))
    return base + gain * (simpleGauss)

def make2DGauss(x,y, mean=0.0, sd=1.0, gain=1.0, base=0.0):
    """
    Return the gaussian distribution for a given set of x-vals

   :Parameters:
        x,y should be x and y indexes  as might be created by numpy.mgrid
        mean: float
            the centre of the distribution - may be a tuple
        sd: float
            the width of the distribution - may be a tuple
        gain: float
            the height of the distribution
        base: float
            an offset added to the result

    """
    
    if numpy.size(sd)==1:
        sd = [sd, sd]
    if numpy.size(mean)==1:
        mean = [mean, mean]
        
    simpleGauss = numpy.exp((-numpy.power(x - mean[0], 2) / (2 * sd[0]**2))-(numpy.power(y - mean[1], 2) / (2 * sd[1]**2)))
    return base + gain * (simpleGauss)

def getRMScontrast(matrix):
    """Returns the RMS contrast (the sample standard deviation) of a array
    """
    RMScontrast = numpy.std(matrix)
    return RMScontrast


def conv2d(smaller, larger):
    """Convolve a pair of 2d numpy matrices.

    Uses fourier transform method, so faster if larger matrix
    has dimensions of size 2**n

    Actually right now the matrices must be the same size (will sort out
    padding issues another day!)
    """
    smallerFFT = fft2(smaller)
    largerFFT = fft2(larger)

    invFFT = ifft2(smallerFFT * largerFFT)
    return invFFT.real


def imfft(X):
    """Perform 2D FFT on an image and center low frequencies
    """
    return fftshift(fft2(X))


def imifft(X):
    """Inverse 2D FFT with decentering
    """
    return numpy.abs(ifft2(ifftshift(X)))


def butter2d_lp(size, cutoff, n=3):
    """Create lowpass 2D Butterworth filter.

       :Parameters:
           size : tuple
               size of the filter
           cutoff : float
               relative cutoff frequency of the filter (0 - 1.0)
           n : int, optional
               order of the filter, the higher n is the sharper
               the transition is.

       :Returns:
           numpy.ndarray
             filter kernel in 2D centered
       """
    if not 0 < cutoff <= 1.0:
        raise ValueError('Cutoff frequency must be between 0 and 1.0')

    if not isinstance(n, int):
        raise ValueError('n must be an integer >= 1')

    rows, cols = size

    x = numpy.linspace(-0.5, 0.5, cols)
    y = numpy.linspace(-0.5, 0.5, rows)

    # An array with every pixel = radius relative to center
    radius = numpy.sqrt((x**2)[numpy.newaxis] + (y**2)[:, numpy.newaxis])

    f = 1 / (1.0 + (radius/cutoff)**(2 * n))   # The filter
    return f


def butter2d_bp(size, cutin, cutoff, n):
    """Bandpass Butterworth filter in two dimensions.

    :Parameters:
        size : tuple
            size of the filter
        cutin : float
            relative cutin  frequency of the filter (0 - 1.0)
        cutoff : float
            relative cutoff frequency of the filter (0 - 1.0)
        n : int, optional
            order of the filter, the higher n is the sharper
            the transition is.

    :Returns:
        numpy.ndarray
          filter kernel in 2D centered

    """

    return butter2d_lp(size, cutoff, n) - butter2d_lp(size, cutin, n)


def butter2d_hp(size, cutoff, n=3):
    """Highpass Butterworth filter in two dimensions.

    :Parameters:
        size : tuple
            size of the filter
        cutoff : float
            relative cutoff frequency of the filter (0 - 1.0)
        n : int, optional
            order of the filter, the higher n is the sharper
            the transition is.

    :Returns:
        numpy.ndarray:
            filter kernel in 2D centered

    """
    return 1.0 - butter2d_lp(size, cutoff, n)


def butter2d_lp_elliptic(size, cutoff_x, cutoff_y, n=3,
                         alpha=0, offset_x=0, offset_y=0):
    """Butterworth lowpass filter of any elliptical shape.

    :Parameters:
        size : tuple
            size of the filter
        cutoff_x, cutoff_y : float, float
            relative cutoff frequency of the filter (0 - 1.0) for x and y axes
        alpha : float, optional
            rotation angle (in radians)
        offset_x, offset_y : float
            offsets for the ellipsoid
        n : int, optional
            order of the filter, the higher n is the sharper
            the transition is.

    :Returns:
        numpy.ndarray:
            filter kernel in 2D centered

    """

    if not (0 < cutoff_x <= 1.0):
        raise ValueError('cutoff_x frequency must be between 0 and 1')
    if not (0 < cutoff_y <= 1.0):
        raise ValueError('cutoff_y frequency must be between 0 and 1')

    rows, cols = size

    # this time we start up with 2D arrays for easy broadcasting
    x = (numpy.linspace(-0.5, 0.5, int(cols)) - offset_x)[numpy.newaxis]
    y = (numpy.linspace(-0.5, 0.5, int(rows)) - offset_y)[:, numpy.newaxis]

    x2 = (x * numpy.cos(alpha) - y * numpy.sin(-alpha))
    y2 = (x * numpy.sin(-alpha) + y * numpy.cos(alpha))

    f = 1 / (1+((x2/(cutoff_x))**2+(y2/(cutoff_y))**2)**n)

    return f

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""Various useful functions for creating filters and textures
(e.g. for PatchStim)
"""

# Part of the PsychoPy library
# Copyright (C) 2002-2018 <NAME> (C) 2019-2021 Open Science Tools Ltd.
# Distributed under the terms of the GNU General Public License (GPL).

from __future__ import absolute_import, division, print_function

from past.utils import old_div
import numpy
from numpy.fft import fft2, ifft2, fftshift, ifftshift
from psychopy import logging
try:
    from PIL import Image
except ImportError:
    import Image


def makeGrating(res,
                ori=0.0,  # in degrees
                cycles=1.0,
                phase=0.0,  # in degrees
                gratType="sin",
                contr=1.0):
    """Make an array containing a luminance grating of the specified params

    :Parameters:
        res: integer
            the size of the resulting matrix on both dimensions (e.g 256)
        ori: float or int (default=0.0)
            the orientation of the grating in degrees
        cycles:float or int (default=1.0)
            the number of grating cycles within the array
        phase: float or int (default=0.0)
            the phase of the grating in degrees (NB this differs to most
            PsychoPy phase arguments which use units of fraction of a cycle)
        gratType: 'sin', 'sqr', 'ramp' or 'sinXsin' (default="sin")
            the type of grating to be 'drawn'
        contr: float (default=1.0)
            contrast of the grating

    :Returns:
        a square numpy array of size resXres

    """
    # to prevent the sinusoid ever being exactly at zero (for sqr wave):
    tiny = 0.0000000000001
    ori *= (old_div(-numpy.pi, 180))
    phase *= (old_div(numpy.pi, 180))
    cyclesTwoPi = cycles * 2.0 * numpy.pi
    xrange, yrange = numpy.mgrid[0.0: cyclesTwoPi: old_div(cyclesTwoPi, res),
                                 0.0: cyclesTwoPi: old_div(cyclesTwoPi, res)]

    sin, cos = numpy.sin, numpy.cos
    if gratType == "none":
        res = 2
        intensity = numpy.ones((res, res), float)
    elif gratType == "sin":
        intensity = contr * sin(xrange * sin(ori) + yrange * cos(ori) + phase)
    elif gratType == "ramp":
        intensity = contr * (xrange * cos(ori) +
                             yrange * sin(ori)) / (2 * numpy.pi)
    elif gratType == "sqr":  # square wave (symmetric duty cycle)
        intensity = numpy.where(sin(xrange * sin(ori) + yrange * cos(ori) +
                                    phase + tiny) >= 0, 1, -1)
    elif gratType == "sinXsin":
        intensity = sin(xrange) * sin(yrange)
    else:
        # might be a filename of an image
        try:
            im = Image.open(gratType)
        except Exception:
            logging.error("couldn't find tex...", gratType)
            return
        # todo: opened it, now what?
    return intensity


def maskMatrix(matrix, shape='circle', radius=1.0, center=(0.0, 0.0)):
    """Make and apply a mask to an input matrix (e.g. a grating)

    :Parameters:
         matrix:  a square numpy array
             array to which the mask should be applied
         shape:  'circle','gauss','ramp' (linear gradient from center)
             shape of the mask
         radius:  float
             scale factor to be applied to the mask (circle with radius of
             [1,1] will extend just to the edge of the matrix). Radius can
             be asymmetric, e.g. [1.0,2.0] will be wider than it is tall.
         center:  2x1 tuple or list (default=[0.0,0.0])
             the centre of the mask in the matrix ([1,1] is top-right
             corner, [-1,-1] is bottom-left)
    """
    # NB makeMask now returns a value -1:1
    alphaMask = makeMask(matrix.shape[0], shape, radius,
                         center=(0.0, 0.0), range=[0, 1])
    return matrix * alphaMask


def makeMask(matrixSize, shape='circle', radius=1.0, center=(0.0, 0.0),
             range=(-1, 1), fringeWidth=0.2):
    """Returns a matrix to be used as an alpha mask (circle,gauss,ramp).

    :Parameters:
            matrixSize: integer
                the size of the resulting matrix on both dimensions (e.g 256)
            shape:  'circle','gauss','ramp' (linear gradient from center),
                'raisedCosine' (the edges are blurred by a raised cosine)
                shape of the mask
            radius:  float
                scale factor to be applied to the mask (circle with radius of
                [1,1] will extend just to the edge of the matrix). Radius can
                asymmetric, e.g. [1.0,2.0] will be wider than it is tall.
            center:  2x1 tuple or list (default=[0.0,0.0])
                the centre of the mask in the matrix ([1,1] is top-right
                corner, [-1,-1] is bottom-left)
            fringeWidth: float (0-1)
                The proportion of the raisedCosine that is being blurred.
            range: 2x1 tuple or list (default=[-1,1])
                The minimum and maximum value in the mask matrix
    """
    rad = makeRadialMatrix(matrixSize, center, radius)
    if shape == 'ramp':
        outArray = 1 - rad
    elif shape == 'circle':
        # outArray=numpy.ones(matrixSize,'f')
        outArray = numpy.where(numpy.greater(rad, 1.0), 0.0, 1.0)
    elif shape == 'gauss':
        outArray = makeGauss(rad, mean=0.0, sd=0.33333)
    elif shape == 'raisedCosine':
        hammingLen = 1000  # This affects the 'granularity' of the raised cos
        fringeProportion = fringeWidth  # This one affects the proportion of
        # the stimulus diameter that is devoted to the raised cosine.

        rad = makeRadialMatrix(matrixSize, center, radius)
        outArray = numpy.zeros_like(rad)
        outArray[numpy.where(rad < 1)] = 1
        raisedCosIdx = numpy.where(
            [numpy.logical_and(rad <= 1, rad >= 1 - fringeProportion)])[1:]

        # Make a raised_cos (half a hamming window):
        raisedCos = numpy.hamming(hammingLen)[:hammingLen//2]
        raisedCos -= numpy.min(raisedCos)
        raisedCos /= numpy.max(raisedCos)

        # Measure the distance from the edge - this is your index into the
        # hamming window:
        dFromEdge = numpy.abs((1 - fringeProportion) - rad[raisedCosIdx])
        dFromEdge /= numpy.max(dFromEdge)
        dFromEdge *= numpy.round(hammingLen/2)

        # This is the indices into the hamming (larger for small distances
        # from the edge!):
        portion_idx = (-1 * dFromEdge).astype(int)

        # Apply the raised cos to this portion:
        outArray[raisedCosIdx] = raisedCos[portion_idx]

        # Sometimes there are some remaining artifacts from this process, get
        # rid of them:
        artifact_idx = numpy.where(
            numpy.logical_and(outArray == 0, rad < 0.99))
        outArray[artifact_idx] = 1
        artifact_idx = numpy.where(
            numpy.logical_and(outArray == 1, rad > 0.99))
        outArray[artifact_idx] = 0

    else:
        raise ValueError('Unknown value for shape argument %s' % shape)
    mag = range[1] - range[0]
    offset = range[0]
    return outArray * mag + offset


def makeRadialMatrix(matrixSize, center=(0.0, 0.0), radius=1.0):
    """Generate a square matrix where each element val is
    its distance from the centre of the matrix

    :Parameters:
        matrixSize: integer
            the size of the resulting matrix on both dimensions (e.g 256)
        radius:  float
            scale factor to be applied to the mask (circle with radius of
            [1,1] will extend just to the edge of the matrix). Radius can
            be asymmetric, e.g. [1.0,2.0] will be wider than it is tall.
        center:  2x1 tuple or list (default=[0.0,0.0])
            the centre of the mask in the matrix ([1,1] is top-right
            corner, [-1,-1] is bottom-left)
    """
    if type(radius) in [int, float]:
        radius = [radius, radius]

    # NB need to add one step length because
    yy, xx = numpy.mgrid[0:matrixSize, 0:matrixSize]
    xx = ((1.0 - 2.0 / matrixSize * xx) + center[0]) / radius[0]
    yy = ((1.0 - 2.0 / matrixSize * yy) + center[1]) / radius[1]
    rad = numpy.sqrt(numpy.power(xx, 2) + numpy.power(yy, 2))
    return rad


def makeGauss(x, mean=0.0, sd=1.0, gain=1.0, base=0.0):
    """
    Return the gaussian distribution for a given set of x-vals

   :Parameters:
        mean: float
            the centre of the distribution
        sd: float
            the width of the distribution
        gain: float
            the height of the distribution
        base: float
            an offset added to the result

    """
    simpleGauss = numpy.exp((-numpy.power(mean - x, 2) / (2 * sd**2)))
    return base + gain * (simpleGauss)

def make2DGauss(x,y, mean=0.0, sd=1.0, gain=1.0, base=0.0):
    """
    Return the gaussian distribution for a given set of x-vals

   :Parameters:
        x,y should be x and y indexes  as might be created by numpy.mgrid
        mean: float
            the centre of the distribution - may be a tuple
        sd: float
            the width of the distribution - may be a tuple
        gain: float
            the height of the distribution
        base: float
            an offset added to the result

    """
    
    if numpy.size(sd)==1:
        sd = [sd, sd]
    if numpy.size(mean)==1:
        mean = [mean, mean]
        
    simpleGauss = numpy.exp((-numpy.power(x - mean[0], 2) / (2 * sd[0]**2))-(numpy.power(y - mean[1], 2) / (2 * sd[1]**2)))
    return base + gain * (simpleGauss)

def getRMScontrast(matrix):
    """Returns the RMS contrast (the sample standard deviation) of a array
    """
    RMScontrast = numpy.std(matrix)
    return RMScontrast


def conv2d(smaller, larger):
    """Convolve a pair of 2d numpy matrices.

    Uses fourier transform method, so faster if larger matrix
    has dimensions of size 2**n

    Actually right now the matrices must be the same size (will sort out
    padding issues another day!)
    """
    smallerFFT = fft2(smaller)
    largerFFT = fft2(larger)

    invFFT = ifft2(smallerFFT * largerFFT)
    return invFFT.real


def imfft(X):
    """Perform 2D FFT on an image and center low frequencies
    """
    return fftshift(fft2(X))


def imifft(X):
    """Inverse 2D FFT with decentering
    """
    return numpy.abs(ifft2(ifftshift(X)))


def butter2d_lp(size, cutoff, n=3):
    """Create lowpass 2D Butterworth filter.

       :Parameters:
           size : tuple
               size of the filter
           cutoff : float
               relative cutoff frequency of the filter (0 - 1.0)
           n : int, optional
               order of the filter, the higher n is the sharper
               the transition is.

       :Returns:
           numpy.ndarray
             filter kernel in 2D centered
       """
    if not 0 < cutoff <= 1.0:
        raise ValueError('Cutoff frequency must be between 0 and 1.0')

    if not isinstance(n, int):
        raise ValueError('n must be an integer >= 1')

    rows, cols = size

    x = numpy.linspace(-0.5, 0.5, cols)
    y = numpy.linspace(-0.5, 0.5, rows)

    # An array with every pixel = radius relative to center
    radius = numpy.sqrt((x**2)[numpy.newaxis] + (y**2)[:, numpy.newaxis])

    f = 1 / (1.0 + (radius/cutoff)**(2 * n))   # The filter
    return f


def butter2d_bp(size, cutin, cutoff, n):
    """Bandpass Butterworth filter in two dimensions.

    :Parameters:
        size : tuple
            size of the filter
        cutin : float
            relative cutin  frequency of the filter (0 - 1.0)
        cutoff : float
            relative cutoff frequency of the filter (0 - 1.0)
        n : int, optional
            order of the filter, the higher n is the sharper
            the transition is.

    :Returns:
        numpy.ndarray
          filter kernel in 2D centered

    """

    return butter2d_lp(size, cutoff, n) - butter2d_lp(size, cutin, n)


def butter2d_hp(size, cutoff, n=3):
    """Highpass Butterworth filter in two dimensions.

    :Parameters:
        size : tuple
            size of the filter
        cutoff : float
            relative cutoff frequency of the filter (0 - 1.0)
        n : int, optional
            order of the filter, the higher n is the sharper
            the transition is.

    :Returns:
        numpy.ndarray:
            filter kernel in 2D centered

    """
    return 1.0 - butter2d_lp(size, cutoff, n)


def butter2d_lp_elliptic(size, cutoff_x, cutoff_y, n=3,
                         alpha=0, offset_x=0, offset_y=0):
    """Butterworth lowpass filter of any elliptical shape.

    :Parameters:
        size : tuple
            size of the filter
        cutoff_x, cutoff_y : float, float
            relative cutoff frequency of the filter (0 - 1.0) for x and y axes
        alpha : float, optional
            rotation angle (in radians)
        offset_x, offset_y : float
            offsets for the ellipsoid
        n : int, optional
            order of the filter, the higher n is the sharper
            the transition is.

    :Returns:
        numpy.ndarray:
            filter kernel in 2D centered

    """

    if not (0 < cutoff_x <= 1.0):
        raise ValueError('cutoff_x frequency must be between 0 and 1')
    if not (0 < cutoff_y <= 1.0):
        raise ValueError('cutoff_y frequency must be between 0 and 1')

    rows, cols = size

    # this time we start up with 2D arrays for easy broadcasting
    x = (numpy.linspace(-0.5, 0.5, int(cols)) - offset_x)[numpy.newaxis]
    y = (numpy.linspace(-0.5, 0.5, int(rows)) - offset_y)[:, numpy.newaxis]

    x2 = (x * numpy.cos(alpha) - y * numpy.sin(-alpha))
    y2 = (x * numpy.sin(-alpha) + y * numpy.cos(alpha))

    f = 1 / (1+((x2/(cutoff_x))**2+(y2/(cutoff_y))**2)**n)

    return f

#!/usr/bin/env python # -*- coding: utf-8 -*- Various useful functions for creating filters and textures (e.g. for PatchStim) # Part of the PsychoPy library # Copyright (C) 2002-2018 <NAME> (C) 2019-2021 Open Science Tools Ltd. # Distributed under the terms of the GNU General Public License (GPL). # in degrees # in degrees Make an array containing a luminance grating of the specified params      :Parameters:         res: integer             the size of the resulting matrix on both dimensions (e.g 256)         ori: float or int (default=0.0)             the orientation of the grating in degrees         cycles:float or int (default=1.0)             the number of grating cycles within the array         phase: float or int (default=0.0)             the phase of the grating in degrees (NB this differs to most             PsychoPy phase arguments which use units of fraction of a cycle)         gratType: 'sin', 'sqr', 'ramp' or 'sinXsin' (default="sin")             the type of grating to be 'drawn'         contr: float (default=1.0)             contrast of the grating      :Returns:         a square numpy array of size resXres # to prevent the sinusoid ever being exactly at zero (for sqr wave): # square wave (symmetric duty cycle) # might be a filename of an image # todo: opened it, now what? Make and apply a mask to an input matrix (e.g. a grating)      :Parameters:          matrix:  a square numpy array              array to which the mask should be applied          shape:  'circle','gauss','ramp' (linear gradient from center)              shape of the mask          radius:  float              scale factor to be applied to the mask (circle with radius of              [1,1] will extend just to the edge of the matrix). Radius can              be asymmetric, e.g. [1.0,2.0] will be wider than it is tall.          center:  2x1 tuple or list (default=[0.0,0.0])              the centre of the mask in the matrix ([1,1] is top-right              corner, [-1,-1] is bottom-left) # NB makeMask now returns a value -1:1 Returns a matrix to be used as an alpha mask (circle,gauss,ramp).      :Parameters:             matrixSize: integer                 the size of the resulting matrix on both dimensions (e.g 256)             shape:  'circle','gauss','ramp' (linear gradient from center),                 'raisedCosine' (the edges are blurred by a raised cosine)                 shape of the mask             radius:  float                 scale factor to be applied to the mask (circle with radius of                 [1,1] will extend just to the edge of the matrix). Radius can                 asymmetric, e.g. [1.0,2.0] will be wider than it is tall.             center:  2x1 tuple or list (default=[0.0,0.0])                 the centre of the mask in the matrix ([1,1] is top-right                 corner, [-1,-1] is bottom-left)             fringeWidth: float (0-1)                 The proportion of the raisedCosine that is being blurred.             range: 2x1 tuple or list (default=[-1,1])                 The minimum and maximum value in the mask matrix # outArray=numpy.ones(matrixSize,'f') # This affects the 'granularity' of the raised cos # This one affects the proportion of # the stimulus diameter that is devoted to the raised cosine. # Make a raised_cos (half a hamming window): # Measure the distance from the edge - this is your index into the # hamming window: # This is the indices into the hamming (larger for small distances # from the edge!): # Apply the raised cos to this portion: # Sometimes there are some remaining artifacts from this process, get # rid of them: Generate a square matrix where each element val is     its distance from the centre of the matrix      :Parameters:         matrixSize: integer             the size of the resulting matrix on both dimensions (e.g 256)         radius:  float             scale factor to be applied to the mask (circle with radius of             [1,1] will extend just to the edge of the matrix). Radius can             be asymmetric, e.g. [1.0,2.0] will be wider than it is tall.         center:  2x1 tuple or list (default=[0.0,0.0])             the centre of the mask in the matrix ([1,1] is top-right             corner, [-1,-1] is bottom-left) # NB need to add one step length because Return the gaussian distribution for a given set of x-vals     :Parameters:         mean: float             the centre of the distribution         sd: float             the width of the distribution         gain: float             the height of the distribution         base: float             an offset added to the result Return the gaussian distribution for a given set of x-vals     :Parameters:         x,y should be x and y indexes  as might be created by numpy.mgrid         mean: float             the centre of the distribution - may be a tuple         sd: float             the width of the distribution - may be a tuple         gain: float             the height of the distribution         base: float             an offset added to the result Returns the RMS contrast (the sample standard deviation) of a array Convolve a pair of 2d numpy matrices.      Uses fourier transform method, so faster if larger matrix     has dimensions of size 2**n      Actually right now the matrices must be the same size (will sort out     padding issues another day!) Perform 2D FFT on an image and center low frequencies Inverse 2D FFT with decentering Create lowpass 2D Butterworth filter.         :Parameters:            size : tuple                size of the filter            cutoff : float                relative cutoff frequency of the filter (0 - 1.0)            n : int, optional                order of the filter, the higher n is the sharper                the transition is.         :Returns:            numpy.ndarray              filter kernel in 2D centered # An array with every pixel = radius relative to center # The filter Bandpass Butterworth filter in two dimensions.      :Parameters:         size : tuple             size of the filter         cutin : float             relative cutin  frequency of the filter (0 - 1.0)         cutoff : float             relative cutoff frequency of the filter (0 - 1.0)         n : int, optional             order of the filter, the higher n is the sharper             the transition is.      :Returns:         numpy.ndarray           filter kernel in 2D centered Highpass Butterworth filter in two dimensions.      :Parameters:         size : tuple             size of the filter         cutoff : float             relative cutoff frequency of the filter (0 - 1.0)         n : int, optional             order of the filter, the higher n is the sharper             the transition is.      :Returns:         numpy.ndarray:             filter kernel in 2D centered Butterworth lowpass filter of any elliptical shape.      :Parameters:         size : tuple             size of the filter         cutoff_x, cutoff_y : float, float             relative cutoff frequency of the filter (0 - 1.0) for x and y axes         alpha : float, optional             rotation angle (in radians)         offset_x, offset_y : float             offsets for the ellipsoid         n : int, optional             order of the filter, the higher n is the sharper             the transition is.      :Returns:         numpy.ndarray:             filter kernel in 2D centered # this time we start up with 2D arrays for easy broadcasting