.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "_auto_examples/shortclips/07_extract_motion_energy.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr__auto_examples_shortclips_07_extract_motion_energy.py: =============================================== Extract motion energy features from the stimuli =============================================== This script describes how to extract motion-energy features from the stimuli. .. Note:: The public data set already contains precomputed motion-energy. Therefore, you do not need to run this script to fit motion-energy models in other part of this tutorial. *Motion-energy features:* Motion-energy features result from filtering a video stimulus with spatio-temporal Gabor filters. A pyramid of filters is used to compute the motion-energy features at multiple spatial and temporal scales. Motion-energy features were introduced in [1]_. The motion-energy extraction is performed by the package `pymoten `_. Check the pymoten `gallery of examples `_ for visualizing motion-energy filters, and for pymoten API usage examples. Running time ------------ Extracting motion energy is a bit longer than the other examples. It typically takes a couple hours to run. .. GENERATED FROM PYTHON SOURCE LINES 27-27 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 29-35 .. code-block:: Python # path of the data directory from voxelwise_tutorials.io import get_data_home directory = get_data_home(dataset="shortclips") print(directory) .. GENERATED FROM PYTHON SOURCE LINES 36-40 Load the stimuli images ----------------------- Here the data is not loaded in memory, we only take a peek at the data shape. .. GENERATED FROM PYTHON SOURCE LINES 40-49 .. code-block:: Python import os import h5py first_file_name = os.path.join(directory, 'stimuli', 'train_00.hdf') print(f"Content of {first_file_name}:") with h5py.File(first_file_name, 'r') as f: for key in f.keys(): print(f[key]) .. GENERATED FROM PYTHON SOURCE LINES 50-58 Compute the luminance --------------------- The motion energy is typically not computed on RGB (color) images, but on the luminance channel of the LAB color space. To avoid loading the entire simulus array in memory, we use batches of data. These batches can be arbitrary, since the luminance is computed independently on each image. .. GENERATED FROM PYTHON SOURCE LINES 58-96 .. code-block:: Python import numpy as np from moten.io import imagearray2luminance from himalaya.progress_bar import bar from voxelwise_tutorials.io import load_hdf5_array def compute_luminance(run_name, size=(96, 96), batch_size=100): stimuli_file = os.path.join(directory, 'stimuli', run_name) # get the number of images in the stimuli file with h5py.File(stimuli_file, 'r') as f: n_images = f['stimuli'].shape[0] # compute the luminance on each batch luminance = np.zeros((n_images, *size)) for start in bar(range(0, n_images, batch_size), title=f'compute_luminance({run_name})'): # load the batch of images batch = slice(start, start + batch_size) images = load_hdf5_array(stimuli_file, key='stimuli', slice=batch) # ``imagearray2luminance`` uses uint8 arrays if images.dtype != 'uint8': images = np.int_(np.clip(images, 0, 1) * 255).astype(np.uint8) # convert RGB images to a single luminance channel luminance[batch] = imagearray2luminance(images, size=size) return luminance luminance_train = np.concatenate( [compute_luminance(f"train_{ii:02d}.hdf") for ii in range(12)]) luminance_test = compute_luminance("test.hdf") .. GENERATED FROM PYTHON SOURCE LINES 97-107 Compute the motion energy ------------------------- This is done with a ``MotionEnergyPyramid`` object of the ``pymoten`` package. The parameters used are the one described in [1]_. Here we use batches corresponding to run lengths. Indeed, motion energy is computed over multiple images, since the filters have a temporal component. Therefore, motion-energy is not independent of other images, and we cannot arbitrarily split the images. .. GENERATED FROM PYTHON SOURCE LINES 107-150 .. code-block:: Python from scipy.signal import decimate from moten.pyramids import MotionEnergyPyramid # fixed experiment settings N_FRAMES_PER_SEC = 15 N_FRAMES_PER_TR = 30 N_TRS_PER_RUN = 300 def compute_motion_energy(luminance, batch_size=N_TRS_PER_RUN * N_FRAMES_PER_TR, noise=0.1): n_frames, height, width = luminance.shape # We create a pyramid instance, with the main motion-energy parameters. pyramid = MotionEnergyPyramid(stimulus_vhsize=(height, width), stimulus_fps=N_FRAMES_PER_SEC, spatial_frequencies=[0, 2, 4, 8, 16, 32]) # We batch images run by run. motion_energy = np.zeros((n_frames, pyramid.nfilters)) for ii, start in enumerate(range(0, n_frames, batch_size)): batch = slice(start, start + batch_size) print("run %d" % ii) # add some noise to deal with constant black areas luminance_batch = luminance[batch].copy() luminance_batch += np.random.randn(*luminance_batch.shape) * noise luminance_batch = np.clip(luminance_batch, 0, 100) motion_energy[batch] = pyramid.project_stimulus(luminance_batch) # decimate to the sampling frequency of fMRI responses motion_energy_decimated = decimate(motion_energy, N_FRAMES_PER_TR, ftype='fir', axis=0) return motion_energy_decimated motion_energy_train = compute_motion_energy(luminance_train) motion_energy_test = compute_motion_energy(luminance_test) .. GENERATED FROM PYTHON SOURCE LINES 151-154 We end this script with saving the features. These features should be approximately equal to the "motion-energy" features already precomputed in the public data set. .. GENERATED FROM PYTHON SOURCE LINES 154-166 .. code-block:: Python from voxelwise_tutorials.io import save_hdf5_dataset features_directory = os.path.join(directory, "features") if not os.path.exists(features_directory): os.makedirs(features_directory) save_hdf5_dataset( os.path.join(features_directory, "motion_energy_recomputed.hdf"), dataset=dict(X_train=motion_energy_train, X_test=motion_energy_test, run_onsets=np.arange(0, 3600, 300))) .. GENERATED FROM PYTHON SOURCE LINES 167-173 References ---------- .. [1] Nishimoto, S., Vu, A. T., Naselaris, T., Benjamini, Y., Yu, B., & Gallant, J. L. (2011). Reconstructing visual experiences from brain activity evoked by natural movies. Current Biology, 21(19), 1641-1646. .. _sphx_glr_download__auto_examples_shortclips_07_extract_motion_energy.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 07_extract_motion_energy.ipynb <07_extract_motion_energy.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 07_extract_motion_energy.py <07_extract_motion_energy.py>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_