Skip to content

Commit e7b745e

Browse files
TomDLTTomDLT
committed
ENH add motion energy extraction script
1 parent b511c1b commit e7b745e

7 files changed

Lines changed: 224 additions & 55 deletions

File tree

‎README.rst‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -44,8 +44,7 @@ Developers can also install the package in editable mode via:
4444
Requirements
4545
------------
4646

47-
The Python package contained in this repository, `voxelwise`, requires the
48-
following dependencies:
47+
The Python package ``voxelwise`` requires the following dependencies:
4948

5049
- `numpy <https://github.com/numpy/numpy>`_
5150
- `scipy <https://github.com/scipy/scipy>`_

‎doc/index.rst‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,8 @@ Welcome to the voxelwise modeling tutorial from the Gallant lab.
88
Tutorials
99
---------
1010
.. toctree::
11-
:maxdepth: 4
11+
:includehidden:
12+
:maxdepth: 2
1213

1314
auto_examples/index
1415

‎tutorials/movies_3T/02_plot_wordnet_model.py‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@
1919
set, comparing the model predictions with the corresponding ground-truth fMRI
2020
responses.
2121
22-
The ridge model uses the package "himalaya", available
22+
The ridge model uses the package ``himalaya``, available
2323
at https://github.com/gallantlab/himalaya.
2424
This package can fit the model either on CPU or on GPU.
2525
"""
@@ -123,8 +123,9 @@
123123
from voxelwise.delayer import Delayer
124124

125125
###############################################################################
126-
# We set the backend to "torch_cuda" to fit the model using GPU.
126+
# We set himalaya's backend to "torch_cuda" to fit the model using GPU.
127127
# The available backends are:
128+
#
128129
# - "numpy" (CPU) (default)
129130
# - "torch" (CPU)
130131
# - "torch_cuda" (GPU)

‎tutorials/movies_3T/03_plot_motion_energy_model.py‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@
2424
set, comparing the model predictions with the corresponding ground-truth fMRI
2525
responses.
2626
27-
The ridge model uses the package "himalaya", available
27+
The ridge model uses the package ``himalaya``, available
2828
at https://github.com/gallantlab/himalaya.
2929
This package can fit the model either on CPU or on GPU.
3030
"""
@@ -129,8 +129,9 @@
129129
from voxelwise.delayer import Delayer
130130

131131
###############################################################################
132-
# We set the backend to "torch_cuda" to fit the model using GPU.
132+
# We set himalaya's backend to "torch_cuda" to fit the model using GPU.
133133
# The available backends are:
134+
#
134135
# - "numpy" (CPU) (default)
135136
# - "torch" (CPU)
136137
# - "torch_cuda" (GPU)

‎tutorials/movies_3T/05_extract_motion_energy.py‎

Lines changed: 157 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,5 +3,162 @@
33
Extract motion energy features from the stimuli
44
===============================================
55
6+
This script describes how to extract motion-energy features from the stimuli.
7+
8+
.. Note::
9+
This public data set already contains precomputed motion-energy. Therefore,
10+
you do not need to run this script to fit motion-energy models in other
11+
part of this tutorial.
12+
13+
Motion-energy
14+
-------------
15+
16+
Motion-energy features result from filtering a video stimulus with
17+
spatio-temporal Gabor filters. A pyramid of filters is used to compute the
18+
motion-energy features at multiple spatial and temporal scales.
19+
Motion-energy features were introduced in [1]_.
20+
21+
The motion-energy extraction is performed by the package ``pymoten``, available
22+
at https://github.com/gallantlab/pymoten.
23+
24+
.. [1] Nishimoto, S., Vu, A. T., Naselaris, T., Benjamini, Y., Yu,
25+
B., & Gallant, J. L. (2011). Reconstructing visual experiences from brain
26+
activity evoked by natural movies. Current Biology, 21(19), 1641-1646.
27+
628
"""
729
# sphinx_gallery_thumbnail_path = "_static/moten.png"
30+
###############################################################################
31+
32+
# path of the data directory
33+
directory = '/data1/tutorials/vim-4/'
34+
35+
###############################################################################
36+
# Load the stimuli images
37+
# -----------------------
38+
#
39+
# Here the data is not loaded in memory, we only take a peak at the data shape.
40+
41+
import h5py
42+
import os
43+
44+
first_file_name = os.path.join(directory, 'stimuli', 'train_00.hdf')
45+
print(f"Content of {first_file_name}:")
46+
with h5py.File(first_file_name, 'r') as f:
47+
for key in f.keys():
48+
print(f[key])
49+
50+
###############################################################################
51+
# Compute the luminance
52+
# ---------------------
53+
#
54+
# The motion energy is typically not computed on RGB (color) images,
55+
# but on the luminance channel of the LAB color space.
56+
# To avoid loading the entire simulus array in memory, we use batches of data.
57+
# These batches can be arbitray, since the luminance is computed independently
58+
# on each image.
59+
60+
import numpy as np
61+
from moten.io import imagearray2luminance
62+
63+
from voxelwise.progress_bar import bar
64+
from voxelwise.io import load_hdf5_array
65+
66+
67+
def compute_luminance(run_name, size=(96, 96)):
68+
69+
stimuli_file = os.path.join(directory, 'stimuli', run_name)
70+
71+
# get the list of batches in the stimuli file
72+
with h5py.File(stimuli_file, 'r') as f:
73+
keys = list(f.keys())
74+
keys.sort() # sort the batches
75+
76+
# compute the luminance on each batch
77+
luminance = []
78+
for key in bar(keys, title=f'compute_luminance({run_name})'):
79+
# load the batch of images
80+
images = load_hdf5_array(stimuli_file, key=key)
81+
82+
# ``imagearray2luminance`` uses uint8 arrays
83+
if images.dtype != 'uint8':
84+
images = np.int_(np.clip(images, 0, 1) * 255).astype(np.uint8)
85+
86+
# convert RGB images to a single luminance channel
87+
luminance.append(imagearray2luminance(images, size=size))
88+
89+
return np.concatenate(luminance)
90+
91+
92+
luminance_train = np.concatenate(
93+
[compute_luminance(f"train_{ii:02d}.hdf") for ii in range(12)])
94+
luminance_test = compute_luminance("test.hdf")
95+
96+
###############################################################################
97+
# Compute the motion energy
98+
# -------------------------
99+
#
100+
# This is done with a ``MotionEnergyPyramid`` object of the ``pymoten``
101+
# package. The parameters used are the one described in [1]_.
102+
#
103+
# Here we use batches corresponding to run lengths. Indeed, motion energy is
104+
# computed over multiple images, since the filters have a temporal component.
105+
# Therefore, motion-energy is not independent of other images, and we cannot
106+
# arbitrarily split the images.
107+
108+
from scipy.signal import decimate
109+
from moten.pyramids import MotionEnergyPyramid
110+
111+
# fixed experiment settings
112+
N_FRAMES_PER_SEC = 15
113+
N_FRAMES_PER_TR = 30
114+
N_TRS_PER_RUN = 300
115+
116+
117+
def compute_motion_energy(luminance,
118+
batch_size=N_TRS_PER_RUN * N_FRAMES_PER_TR,
119+
noise=0.1):
120+
121+
n_frames, height, width = luminance.shape
122+
123+
# We create a pyramid instance, with the main motion-energy parameters.
124+
pyramid = MotionEnergyPyramid(stimulus_vhsize=(height, width),
125+
stimulus_fps=N_FRAMES_PER_SEC,
126+
spatial_frequencies=[0, 2, 4, 8, 16, 32])
127+
128+
# We batch images run by run.
129+
motion_energy = np.zeros((n_frames, pyramid.nfilters))
130+
for ii, start in enumerate(range(0, n_frames, batch_size)):
131+
batch = slice(start, start + batch_size)
132+
print("run %d" % ii)
133+
134+
# add some noise to deal with constant black areas
135+
luminance_batch = luminance[batch].copy()
136+
luminance_batch += np.random.randn(*luminance_batch.shape) * noise
137+
luminance_batch[luminance_batch < 0] = 0
138+
luminance_batch[luminance_batch > 100] = 100
139+
140+
motion_energy[batch] = pyramid.project_stimulus(luminance_batch)
141+
142+
# decimate to the sampling frequency of fMRI responses
143+
motion_energy_decimated = decimate(motion_energy, N_FRAMES_PER_TR,
144+
ftype='fir', axis=0)
145+
return motion_energy_decimated
146+
147+
148+
motion_energy_train = compute_motion_energy(luminance_train)
149+
motion_energy_test = compute_motion_energy(luminance_test)
150+
151+
###############################################################################
152+
# We end this script with saving the features. These features should be
153+
# approximately equal to the "motion_energy" features already precomputed in
154+
# the public data set.
155+
156+
from voxelwise.io import save_hdf5_dataset
157+
158+
features_directory = os.path.join(directory, "features")
159+
if not os.path.exists(features_directory):
160+
os.makedirs(features_directory)
161+
162+
save_hdf5_dataset(
163+
os.path.join(features_directory, "motion_energy_recomputed.hdf"),
164+
dataset=dict(X_train=motion_energy_train, X_test=motion_energy_test))

‎tutorials/movies_4T/01_extract_motion_energy.py‎

Lines changed: 31 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -5,40 +5,49 @@
55
66
This script describes how to extract motion-energy features from the stimuli.
77
8+
Motion-energy
9+
-------------
10+
811
Motion-energy features result from filtering a video stimulus with
912
spatio-temporal Gabor filters. A pyramid of filters is used to compute the
1013
motion-energy features at multiple spatial and temporal scales.
14+
Motion-energy features were introduced in [1]_.
1115
12-
The motion-energy extraction is performed by the package "pymoten", available
16+
The motion-energy extraction is performed by the package ``pymoten``, available
1317
at https://github.com/gallantlab/pymoten.
1418
19+
.. [1] Nishimoto, S., Vu, A. T., Naselaris, T., Benjamini, Y., Yu,
20+
B., & Gallant, J. L. (2011). Reconstructing visual experiences from brain
21+
activity evoked by natural movies. Current Biology, 21(19), 1641-1646.
22+
1523
"""
1624
# sphinx_gallery_thumbnail_path = "_static/moten.png"
1725
###############################################################################
18-
26+
# Load the stimuli images
27+
# -----------------------
28+
#
1929
# We downloaded the files in the previous script, and here we update the path
2030
# variable to link to the directory containing the data.
2131

2232
directory = '/data1/tutorials/vim-2/'
2333

2434
###############################################################################
25-
# Then, we preload the stimuli.
26-
#
2735
# Here the data is not loaded in memory, we only take a peak at the data shape.
2836

2937
import h5py
30-
import os.path as op
38+
import os
3139

32-
with h5py.File(op.join(directory, 'Stimuli.mat'), 'r') as f:
40+
with h5py.File(os.path.join(directory, 'Stimuli.mat'), 'r') as f:
3341
print(f.keys()) # Show all variables
3442

3543
for key in f.keys():
3644
print(f[key])
3745

3846
###############################################################################
39-
# Then, we compute the luminance of the stimulus images.
47+
# Compute the luminance
48+
# ---------------------
4049
#
41-
# Indeed, the motion energy is typically not computed on RGB (color) images,
50+
# The motion energy is typically not computed on RGB (color) images,
4251
# but on the luminance channel of the LAB color space.
4352
# To avoid loading the entire simulus array in memory, we use batches of data.
4453
# These batches can be arbitray, since the luminance is computed independently
@@ -52,7 +61,7 @@
5261

5362
def compute_luminance(train_or_test, batch_size=1024):
5463

55-
with h5py.File(op.join(directory, 'Stimuli.mat'), 'r') as f:
64+
with h5py.File(os.path.join(directory, 'Stimuli.mat'), 'r') as f:
5665

5766
if train_or_test == 'train':
5867
data = f['st']
@@ -85,10 +94,11 @@ def compute_luminance(train_or_test, batch_size=1024):
8594
luminance_test = compute_luminance("test")
8695

8796
###############################################################################
88-
# Finally, we compute the motion energy features.
97+
# Compute the motion energy
98+
# -------------------------
8999
#
90-
# This is done with a `MotionEnergyPyramid` object of the `pymoten` package.
91-
# The parameters used are the one described in [Nishimoto et al. 2011].
100+
# This is done with a ``MotionEnergyPyramid`` object of the ``pymoten``
101+
# package. The parameters used are the one described in [1]_.
92102
#
93103
# Here we use batches corresponding to run lengths. Indeed, motion energy is
94104
# computed over multiple images, since the filters have a temporal component.
@@ -139,14 +149,15 @@ def compute_motion_energy(luminance,
139149
motion_energy_test = compute_motion_energy(luminance_test)
140150

141151
###############################################################################
142-
# We end with saving the features.
152+
# We end this script with saving the features, to use them in voxelwise
153+
# modeling in the following example.
143154

144-
import os
155+
from voxelwise.io import save_hdf5_dataset
145156

146-
features_directory = op.join(directory, "features")
147-
if not op.exists(features_directory):
157+
features_directory = os.path.join(directory, "features")
158+
if not os.path.exists(features_directory):
148159
os.makedirs(features_directory)
149-
np.save(op.join(features_directory, "motion_energy_train.npy"),
150-
motion_energy_train)
151-
np.save(op.join(features_directory, "motion_energy_test.npy"),
152-
motion_energy_test)
160+
161+
save_hdf5_dataset(
162+
os.path.join(features_directory, "motion_energy.hdf"),
163+
dataset=dict(X_train=motion_energy_train, X_test=motion_energy_test))

0 commit comments

Comments
 (0)