 # Source code for psychopy.visual.filters

```#!/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 Jonathan Peirce (C) 2019-2022 Open Science Tools Ltd.

import numpy
from numpy.fft import fft2, ifft2, fftshift, ifftshift

[docs]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 *= -numpy.pi / 180.
phase *= numpy.pi / 180.
cyclesTwoPi = cycles * 2.0 * numpy.pi
xrange, yrange = numpy.mgrid[
0.0:cyclesTwoPi:(cyclesTwoPi / res),
0.0:cyclesTwoPi:(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?
raise ValueError("Invalid value for parameter `gratType`.")

return intensity

"""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)
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
center=(0.0, 0.0), range=[0, 1])

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)
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
"""
if shape == 'ramp':
elif shape == 'circle':
# outArray=numpy.ones(matrixSize,'f')
outArray = numpy.where(numpy.greater(rad, 1.0), 0.0, 1.0)
elif shape == 'gauss':
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.

raisedCosIdx = numpy.where(

# 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 - range
offset = range
return outArray * mag + offset

"""Generate a square matrix where each element values is its distance from
the centre of the matrix.

Parameters
----------
matrixSize : int
Matrix size. Corresponds to the number of elements along each dimension.
Must be >0.
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)

Returns
-------
ndarray
Square matrix populated with distance values and
`size == (matrixSize, matrixSize)`.

"""

try:
matrixSize = int(matrixSize)
except ValueError:
raise TypeError('parameter `matrixSize` must be a numeric type')

if matrixSize <= 1:
raise ValueError(
'parameter `matrixSize` must be positive and greater than 1, got: {}'.format(
matrixSize))

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

[docs]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

[docs]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, 2) / (2 * sd**2))-(numpy.power(y - mean, 2) / (2 * sd**2)))
return base + gain * simpleGauss

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

[docs]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
"""
smallerFFT = fft2(smaller)
largerFFT = fft2(larger)

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

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

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

[docs]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

[docs]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)

[docs]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)

[docs]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
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
```