In [1]:
import sys
sys.path.append('../../pyrtools/')
import pyrtools as pt
import numpy as np
sys.path.append('..')
import sfp
%matplotlib inline
import seaborn as sns
import pandas as pd
from scipy import signal
import fractions
import matplotlib.pyplot as plt
import pickle
%load_ext autoreload
%autoreload 2
The following demonstrates some of the stimuli that can be made with the log_polar_grating command. The key characteristic of these stimuli is that the frequency (in cycles per pixel) increases as you approach the center of the image. This encodes our hypothesis that the preferred spatial frequency of neurons (and thus, voxels) whose receptive fields are located at the fovea (here, the center of the image) will have a higher preferred spatial frequency than those whose receptive fields are located in the periphery.
However, with these parameters, we notice that there is strange-looking aliasing happening in the center of the image.
In [4]:
ims = []
freq=2
phase = 0
ampl=1
shape=64
origin = None
ims.append(sfp.stimuli.log_polar_grating(shape, freq,phi=phase,ampl=ampl, origin=origin))
ims.append(sfp.stimuli.log_polar_grating(shape, w_a=freq,phi=phase,ampl=ampl))
ims.append(sfp.stimuli.log_polar_grating(shape, freq/2,freq,phi=phase,ampl=ampl))
ims.append(sfp.stimuli.log_polar_grating(shape,-freq,freq,phi=phase,ampl=ampl))
ims.append(sfp.stimuli.log_polar_grating(shape,freq/2,freq,phi=phase,ampl=ampl))
ims.append(sfp.stimuli.log_polar_grating(shape,freq,freq/4,phi=phase,ampl=ampl))
ims.append(ims[0]+ims[1])
# ims.append(sfp.stimuli._create_better_sampled_grating(shape, freq, phi=phase, ampl=ampl, orig_origin=origin, check_scale_factor=3))
# ims.append(sfp.stimuli._create_better_sampled_grating(shape, w_a=freq, phi=phase, ampl=ampl, orig_origin=origin, check_scale_factor=3))
# ims.append(sfp.stimuli._create_better_sampled_grating(shape, freq, phi=phase, ampl=ampl, check_scale_factor=50, orig_origin=origin))
# ims.append(sfp.stimuli._create_better_sampled_grating(shape, freq, phi=phase, ampl=ampl, check_scale_factor=200, orig_origin=origin))
# Example stimuli
pt.imshow(ims, col_wrap=min(len(ims), 4),zoom=5);
In order to investigate this a little more, we create compare the second, radial stimulus with a more heavily-sampled version.
The first plot shows the central slice of the 64 x 64 image in orange and the central slice of the grating 11 times larger in blue. We clearly see some very high frequency signal that is being missed.
In [14]:
stim, over_sampled_stim = sfp.stimuli.check_aliasing(64, 0, w_a=6, check_scale_factor=11)
This second plot shows the two images, and again we can see the aliasing at the center
In [21]:
pt.imshow([stim, over_sampled_stim], zoom=6/11);
To deal with this issue, we create a mask that we'll lay at the center of the image, completely masking out the aliased portion of the stimuli and then fading gradually to invisible as we move away from it.
This first plot shows the same slice as above for the un-masked, fade-masked, and binary-masked stimuli.
In [24]:
tmps = sfp.stimuli.check_aliasing_with_mask(64, 0, w_a=6, )
This here shows the 64 x 64 image with no mask, the faded mask, and the binary mask. We can see that the masking hides the aliased portion of the image.
In [31]:
pt.imshow([tmps[0], tmps[1] * tmps[0], tmps[2] * tmps[0]], col_wrap=3, zoom=5);
Now that we know the stimuli we want to create and are convinced that masking takes care of their aliasing, we can create the set of stimuli for experiment. We do this using the gen_stim_set function, which requires that the size, origin, and number_of_fade_pixels parameters be constant across all stimuli, but allows the others to vary. We create all the stimuli and then apply the largest mask to them, so that all stimuli have the same mask and none show any aliasing. The function returns both the unmasked (stim) and masked (mstim) versions of the stimuli, but we only want the masked versions (the unmasked are only returned so they can be double-checked).
The following two plots show an example set of (smaller) stimuli and one full-sized one.
To create the exact set of stimuli wanted for the experiment, use the stimuli.main function (but note that this can take some time and will save the outputs; it can also be run on the command-line).
In [37]:
mstim, stim, mask = sfp.stimuli.gen_log_polar_stim_set(32, freqs_ra=[(4,0),(0,4),(4,4)], phi=np.array(range(8))/8.*2*np.pi, )
pt.imshow(mstim, col_wrap=6, zoom=5);
In [40]:
mstim, _, _ = sfp.stimuli.gen_log_polar_stim_set(1080, [(0, 128)], bytescale=True)# phi=np.array(range(10))/10.*2*np.pi, )
pt.imshow(mstim, col_wrap=4, );
For the purposes of creating our stimuli, we're interested in a couple of numbers that aren't apparent from the plots above:
for the spatial frequencies, we're interested in this in units of cycles per visual degree, so we create the full-size stimulus (1080 x 1080 pixels) and input the diameter of the stimulus in visual degrees (28).
In [60]:
# Note that this takes a while to run
mask_df, sf_df = sfp.stimuli.check_stim_properties(1080, None, 24, w_r=range(0,181,5), w_a=range(0,181,5))
We see that the mask grows with the frequencies. It gets pretty big.
Note that the mask size is only dependent on the frequency content the image, not the size of the image. So if you have $\omega_a=500$, your mask will always be 186 pixels wide, so if your image is only 150 pixels per side, then your entire image will be blank.
In [70]:
sfp.stimuli.plot_stim_properties(mask_df, size=8, data_label='mask_radius_pix', title_text='Mask radius in degrees')
This next plot isn't as interesting, but it's a sanity check to make sure the mask is working: basically, we want to see that the mask is working and so nowhere do we get frequencies exceeding .5 cycles per pixel. This shows that that's the case: no matter how high the frequencies go in either direction, they always stay below .5
In [71]:
sfp.stimuli.plot_stim_properties(mask_df, data_label='cpp_masked_max', title_text="Max masked frequency in cpp", size=8)
What we want to do is ensure that our mask has a radius of about one degree of visual angle, so we plot the distance from the origin in frequency space (that is, $\sqrt{\omega_\alpha^2+\omega_r^2}$) against the radius of the mask, in degrees. Because of how we constructed these stimuli, the local spatial frequency will only depend on this distance (i.e., it's the same in all directions: the stimuli with $\omega_\alpha=0, \omega_r=100$ has the same local spatial frequency everywhere as the one with $\omega_\alpha=100, \omega_r=0$) and the radius of the mask depends on this local spatial frequency (it will mask everywhere the local spatial frequency exceeds our Nyquist frequency limit, which I've set to a slightly conservative .475), this plot will be a single straight line (if you add the hue argument that's currently commented out, you'll see that they all lie on top of each other).
What we want to do is find the distance in frequency space where we cross the mask_radius_deg=1 line, which I show below the figure
In [73]:
g = sns.FacetGrid(mask_df, size=7,)# hue='w_r')
g.map(plt.scatter, 'freq_distance', 'mask_radius_deg')
#sns.plt.scatter(mask_df['freq_distance'], mask_df['mask_radius_deg'])
g.add_legend()
g.ax.plot(g.ax.get_xlim(), [1, 1], 'k--')
g.ax.set_xlim((0, 300))
_=g.ax.set_ylim((0, 2))
In [77]:
mask_df.iloc[abs(mask_df.mask_radius_deg-1).argmin()]
Out[77]:
Thus we want to go out to a frequency distance of about 135. Since we're $\log_2$ spacing our stimuli, we want to know the following and see that our stimuli should lie at most about $2^7$ from the origin of our frequency space.
In [75]:
np.log2(134.164079)
Out[75]:
We also want to create stimuli whose local spatial frequency is the same everywhere, i.e., normal 2d sine waves. To compare them effectively to our stimuli, we construct the spatial frequency dataframe.
Because only the distance from the origin in frequency space matters (as discussed above), we set $\omega_r$ to 0 and only look at different values of $\omega_a$, for simplicity. We also use the actual values we will use in this experiment: $\omega_a=[2^{2.5}, 2^3, 2^{3.5}, 2^4, 2^{4.5}, 2^5, 2^{5.5}, 2^6, 2^{6.5}, 2^7]$. We use these values rounded to the nearest integer, because they need to be integers in order to not have weird discontinuities
In order to determine what spatial frequencies are constant grating stimuli should have, we look at the local period (which is equivalent to looking at the local spatial frequency, since they're just reciprocals of each other, and more interpretable) at eccentricity 4.5 degrees. If we construct a constant grating with those periods, we guarantee that at one eccentricity we have the same local spatial frequency in both sets of stimuli. We choose 4.5 because then there are roughly equal numbers of larger and smaller log polar gratings.
In [79]:
_, sf_df = sfp.stimuli.check_stim_properties(1080, None, 24, 0, np.round([2**i for i in np.arange(2.5, 7.5, .5)]))
In [82]:
plot_func = [plt.plot, plt.semilogy]
plot_kwargs = [{}, {'basey': 2}]
plot_idx = 1
g = sns.FacetGrid(sf_df, hue='w_a', size=7)
g.map(plot_func[plot_idx], 'eccen', 'local_period_ppc', **plot_kwargs[plot_idx])
g.add_legend()
constant_grads = sf_df[sf_df.eccen==4.5].local_period_ppc.values
for val in constant_grads:
plot_func[plot_idx]([1.5, 11.5], [val, val], 'k--', **plot_kwargs[plot_idx])
plot_func[plot_idx]([4.5, 4.5], [2**1.5, 2**9], 'k--', **plot_kwargs[plot_idx])
Out[82]:
Unfortunately, these numbers are ugly, but they're the ones we want to use.
In [84]:
print(1./constant_grads)
print(np.log2(1./constant_grads))
In [85]:
stim = np.load("../data/stimuli/unshuffled.npy")
df = pd.read_csv("../data/stimuli/unshuffled_stim_description.csv")
constant_stim = np.load("../data/stimuli/constant_unshuffled.npy")
constant_df = pd.read_csv("../data/stimuli/constant_unshuffled_stim_description.csv")
In [86]:
pt.imshow([stim[0], stim[100], stim[354]]);
In [87]:
pt.imshow([constant_stim[0], constant_stim[100], constant_stim[300]]);
In [88]:
pt.imshow([i for i in stim[:10]], col_wrap=4);
In [90]:
pt.imshow([i for i in stim[300:304]], col_wrap=4);
Let's save some of these stimuli
In [91]:
w_y = constant_df[constant_df.w_x==0].w_y.unique()
stim_idx = constant_df[(constant_df.w_x==0)&(constant_df.w_y.isin(w_y))&(constant_df.phi==0)].index
pt.imshow([i for i in constant_stim[stim_idx]], col_wrap=4);
In [97]:
fig = pt.imshow(stim[:8,:,:],zoom=180/1080,)
fig.savefig('stim_phases.svg')
At a certain point during this experiment, we realized that the scanner projector, like an display system, will lose contrast for higher spatial frequencies. Since it has a pointspread function, this will act as a small amount of blur, which will reduce contrast as the spatial frequency increases. In a side project, we measured the Modulation Transfer Function of the projector for different presentation frequencies. We now load in that MTF (as a spline object) and create stimuli that invert it. Because the values of the stimuli must lie between 0 and 255 in order to be displayed, we must decrease the amplitude for low frequencies (because we cannot increase the amplitude for high frequencies) in order to present all frequencies at the same contrast
In [3]:
with open('/mnt/prince_scratch/spatial_frequency_preferences/stimuli/mtf_func.pkl', 'rb') as f:
mtf = pickle.load(f)
In [9]:
sfs = np.linspace(0, 1)
plt.semilogx(sfs, mtf(sfs),basex=2)
plt.xlabel('spatial frequency (cpp)')
plt.ylabel('contrast')
plt.title('MTF of the display')
Out[9]:
In [75]:
mask = []
stim = []
masked_stim = []
rescaled_stim = []
for f in [2**5, 2**7]:
stim.append(sfp.stimuli.log_polar_grating(1080, f, 0))
_, _, sf_map, _ = sfp.stimuli.create_sf_maps_cpp(1080, w_r=f, w_a=0)
# this is the max frequency we generate, so we use it here
mask.append(sfp.stimuli.create_antialiasing_mask(1080, 128, 0)[1])
mask.append(sfp.stimuli.create_outer_mask(1080, None)[1])
anti_mtf_map = 1./mtf(sf_map)
rescaled_stim.append(anti_mtf_map * stim[-1])
mask = np.logical_and.reduce(mask)
mask = sfp.stimuli._fade_mask(mask, 3, 3, None)
masked_stim = np.array([mask*s for s in stim])
rescaled_stim = np.array([mask*s for s in rescaled_stim])
masked_stim = sfp.utils.bytescale(masked_stim)
rescaled_stim = sfp.utils.bytescale(rescaled_stim)
In [77]:
pt.imshow([mask*sf_map, mtf(mask*sf_map), 1./mtf(mask*sf_map)], zoom=.5);
Notice that the lower frequency grating no longer uses the full dynamic range [0, 255]. This is how it should work.
In [95]:
pt.imshow(np.concatenate([masked_stim, rescaled_stim]), 'auto1', zoom=.5, col_wrap=2);
In [94]:
fig, axes = plt.subplots(1,2,figsize=(30,10))
for i in range(2):
axes[i].plot(masked_stim[i,:,540])
axes[i].plot(rescaled_stim[i,:,540], '--')
We do this for all stimuli, which we'll now load and and look at several mid-slices
In [141]:
rescaled_stim = np.load('/mnt/winawerlab/Projects/spatial_frequency_preferences/BIDS/stimuli/task-sfprescaled_stimuli.npy')
stim = np.load('/mnt/winawerlab/Projects/spatial_frequency_preferences/BIDS/stimuli/task-sfp_stimuli.npy')
stim_df = pd.read_csv('/mnt/winawerlab/Projects/spatial_frequency_preferences/BIDS/stimuli/task-sfp_stim_description.csv')
stim_df = stim_df.drop_duplicates('class_idx')
The following plots are a little strange, since this isn't always the best way to slice through one of these gratings. You'll also notice that, for relatively low frequencies, there's not a huge difference in the amount they're re-scaled. This is because of the shape of the MTF. You don't start to really lose contrast until the frequency gets higher, so the relatively low frequencies are all presented at roughly the same contrast and, thus, their amplitudes are rescaled by roughly the same amount.
In [177]:
fig, axes = plt.subplots(6, 8, True, True, figsize=(40, 30))
for i, (ax, s, rs) in enumerate(zip(axes.flatten(), stim[:48*8:8], rescaled_stim[:48*8:8])):
ax.plot(s[:,540])
ax.plot(rs[:,540], '--')
ax.set_title('class_idx %02d' % i)
plt.savefig('rescaled_stim.svg')