Note

Go to the end to download the full example code. or to run this example in your browser via Binder

Reconstruction of visual stimuli from Miyawaki et al. 2008¶

This example reproduces the experiment presented in Miyawaki et al.[1].

It reconstructs 10x10 binary images from functional MRI data. Random images are used as training set and structured images are used for reconstruction.

The code is a bit elaborate as the example uses, as the original article, a multiscale prediction on the images seen by the subject.

For an encoding approach for the same dataset, see also Encoding models for visual stimuli from Miyawaki et al. 2008

Note

If you are using Nilearn with a version older than 0.9.0, then you should either upgrade your version or import maskers from the input_data module instead of the maskers module.

That is, you should manually replace in the following example all occurrences of:

from nilearn.maskers import NiftiMasker

with:

from nilearn.input_data import NiftiMasker

import sys
import time

First we load the Miyawaki dataset¶

from nilearn import datasets
from nilearn.plotting import show

sys.stderr.write("Fetching dataset...")
t0 = time.time()

miyawaki_dataset = datasets.fetch_miyawaki2008()

# print basic information on the dataset
print(
    "First functional nifti image (4D) is located "
    f"at: {miyawaki_dataset.func[0]}"
)

X_random_filenames = miyawaki_dataset.func[12:]
X_figure_filenames = miyawaki_dataset.func[:12]
y_random_filenames = miyawaki_dataset.label[12:]
y_figure_filenames = miyawaki_dataset.label[:12]
y_shape = (10, 10)

sys.stderr.write(f" Done ({time.time() - t0:.2f}s).\n")

Fetching dataset...[get_dataset_dir] Dataset found in /home/runner/nilearn_data/miyawaki2008
First functional nifti image (4D) is located at: /home/runner/nilearn_data/miyawaki2008/func/data_figure_run01.nii.gz
 Done (0.18s).

Then we prepare and mask the data¶

import numpy as np

from nilearn.maskers import MultiNiftiMasker

sys.stderr.write("Preprocessing data...")
t0 = time.time()

# Load and mask fMRI data
masker = MultiNiftiMasker(
    mask_img=miyawaki_dataset.mask, detrend=True, standardize=False, n_jobs=2
)
masker.fit()
X_train = masker.transform(X_random_filenames)
X_test = masker.transform(X_figure_filenames)

y_train = [
    np.reshape(
        np.loadtxt(y, dtype=int, delimiter=","), (-1, *y_shape), order="F"
    )
    for y in y_random_filenames
]
y_test = [
    np.reshape(
        np.loadtxt(y, dtype=int, delimiter=","), (-1, *y_shape), order="F"
    )
    for y in y_figure_filenames
]
X_train = np.vstack([x[2:] for x in X_train])
y_train = np.vstack([y[:-2] for y in y_train]).astype(float)
X_test = np.vstack([x[2:] for x in X_test])
y_test = np.vstack([y[:-2] for y in y_test]).astype(float)

n_features = X_train.shape[1]


def flatten(list_of_2d_array):
    flattened = [array.ravel() for array in list_of_2d_array]
    return flattened


# Build the design matrix for multiscale computation
# Matrix is squared, y_rows == y_cols
y_cols = y_shape[1]

# Original data
design_matrix = np.eye(100)

# Example of matrix used for multiscale (sum pixels vertically)
#
# 0.5 *
#
# 1 1 0 0 0 0 0 0 0 0
# 0 1 1 0 0 0 0 0 0 0
# 0 0 1 1 0 0 0 0 0 0
# 0 0 0 1 1 0 0 0 0 0
# 0 0 0 0 1 1 0 0 0 0
# 0 0 0 0 0 1 1 0 0 0
# 0 0 0 0 0 0 1 1 0 0
# 0 0 0 0 0 0 0 1 1 0
# 0 0 0 0 0 0 0 0 1 1

height_tf = (np.eye(y_cols) + np.eye(y_cols, k=1))[: y_cols - 1] * 0.5
width_tf = height_tf.T

yt_tall = [np.dot(height_tf, m) for m in y_train]
yt_large = [np.dot(m, width_tf) for m in y_train]
yt_big = [np.dot(height_tf, np.dot(m, width_tf)) for m in y_train]

# Add it to the training set
y_train = [
    np.r_[y.ravel(), t.ravel(), l.ravel(), b.ravel()]
    for y, t, l, b in zip(y_train, yt_tall, yt_large, yt_big)
]

y_test = np.asarray(flatten(y_test))
y_train = np.asarray(y_train)

# Remove rest period
X_train = X_train[y_train[:, 0] != -1]
y_train = y_train[y_train[:, 0] != -1]
X_test = X_test[y_test[:, 0] != -1]
y_test = y_test[y_test[:, 0] != -1]

sys.stderr.write(f" Done ({time.time() - t0:.2f}s).\n")

Preprocessing data...[MultiNiftiMasker.fit] Resampling mask
 Done (13.34s).

We define our prediction function¶

sys.stderr.write("Training classifiers... \r")
t0 = time.time()

from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.linear_model import OrthogonalMatchingPursuit
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

# Create as many OrthogonalMatchingPursuit as voxels to predict
clfs = []
n_clfs = y_train.shape[1]
for i in range(y_train.shape[1]):
    sys.stderr.write(
        f"Training classifiers {int(i + 1):03}/{int(n_clfs)}... \r"
    )

    clf = Pipeline(
        [
            ("selection", SelectKBest(f_classif, k=500)),
            ("scl", StandardScaler()),
            ("clf", OrthogonalMatchingPursuit(n_nonzero_coefs=10)),
        ]
    )
    clf.fit(X_train, y_train[:, i])
    clfs.append(clf)

sys.stderr.write(
    f"Training classifiers {n_clfs:03d}/{n_clfs:d}... "
    f"Done ({(time.time() - t0):.2f}s).\n"
)

Training classifiers...
Training classifiers 001/361...
Training classifiers 002/361...
Training classifiers 003/361...
Training classifiers 004/361...
Training classifiers 005/361...
Training classifiers 006/361...
Training classifiers 007/361...
Training classifiers 008/361...
Training classifiers 009/361...
Training classifiers 010/361...
Training classifiers 011/361...
Training classifiers 012/361...
Training classifiers 013/361...
Training classifiers 014/361...
Training classifiers 015/361...
Training classifiers 016/361...
Training classifiers 017/361...
Training classifiers 018/361...
Training classifiers 019/361...
Training classifiers 020/361...
Training classifiers 021/361...
Training classifiers 022/361...
Training classifiers 023/361...
Training classifiers 024/361...
Training classifiers 025/361...
Training classifiers 026/361...
Training classifiers 027/361...
Training classifiers 028/361...
Training classifiers 029/361...
Training classifiers 030/361...
Training classifiers 031/361...
Training classifiers 032/361...
Training classifiers 033/361...
Training classifiers 034/361...
Training classifiers 035/361...
Training classifiers 036/361...
Training classifiers 037/361...
Training classifiers 038/361...
Training classifiers 039/361...
Training classifiers 040/361...
Training classifiers 041/361...
Training classifiers 042/361...
Training classifiers 043/361...
Training classifiers 044/361...
Training classifiers 045/361...
Training classifiers 046/361...
Training classifiers 047/361...
Training classifiers 048/361...
Training classifiers 049/361...
Training classifiers 050/361...
Training classifiers 051/361...
Training classifiers 052/361...
Training classifiers 053/361...
Training classifiers 054/361...
Training classifiers 055/361...
Training classifiers 056/361...
Training classifiers 057/361...
Training classifiers 058/361...
Training classifiers 059/361...
Training classifiers 060/361...
Training classifiers 061/361...
Training classifiers 062/361...
Training classifiers 063/361...
Training classifiers 064/361...
Training classifiers 065/361...
Training classifiers 066/361...
Training classifiers 067/361...
Training classifiers 068/361...
Training classifiers 069/361...
Training classifiers 070/361...
Training classifiers 071/361...
Training classifiers 072/361...
Training classifiers 073/361...
Training classifiers 074/361...
Training classifiers 075/361...
Training classifiers 076/361...
Training classifiers 077/361...
Training classifiers 078/361...
Training classifiers 079/361...
Training classifiers 080/361...
Training classifiers 081/361...
Training classifiers 082/361...
Training classifiers 083/361...
Training classifiers 084/361...
Training classifiers 085/361...
Training classifiers 086/361...
Training classifiers 087/361...
Training classifiers 088/361...
Training classifiers 089/361...
Training classifiers 090/361...
Training classifiers 091/361...
Training classifiers 092/361...
Training classifiers 093/361...
Training classifiers 094/361...
Training classifiers 095/361...
Training classifiers 096/361...
Training classifiers 097/361...
Training classifiers 098/361...
Training classifiers 099/361...
Training classifiers 100/361...
Training classifiers 101/361...
Training classifiers 102/361...
Training classifiers 103/361...
Training classifiers 104/361...
Training classifiers 105/361...
Training classifiers 106/361...
Training classifiers 107/361...
Training classifiers 108/361...
Training classifiers 109/361...
Training classifiers 110/361...
Training classifiers 111/361...
Training classifiers 112/361...
Training classifiers 113/361...
Training classifiers 114/361...
Training classifiers 115/361...
Training classifiers 116/361...
Training classifiers 117/361...
Training classifiers 118/361...
Training classifiers 119/361...
Training classifiers 120/361...
Training classifiers 121/361...
Training classifiers 122/361...
Training classifiers 123/361...
Training classifiers 124/361...
Training classifiers 125/361...
Training classifiers 126/361...
Training classifiers 127/361...
Training classifiers 128/361...
Training classifiers 129/361...
Training classifiers 130/361...
Training classifiers 131/361...
Training classifiers 132/361...
Training classifiers 133/361...
Training classifiers 134/361...
Training classifiers 135/361...
Training classifiers 136/361...
Training classifiers 137/361...
Training classifiers 138/361...
Training classifiers 139/361...
Training classifiers 140/361...
Training classifiers 141/361...
Training classifiers 142/361...
Training classifiers 143/361...
Training classifiers 144/361...
Training classifiers 145/361...
Training classifiers 146/361...
Training classifiers 147/361...
Training classifiers 148/361...
Training classifiers 149/361...
Training classifiers 150/361...
Training classifiers 151/361...
Training classifiers 152/361...
Training classifiers 153/361...
Training classifiers 154/361...
Training classifiers 155/361...
Training classifiers 156/361...
Training classifiers 157/361...
Training classifiers 158/361...
Training classifiers 159/361...
Training classifiers 160/361...
Training classifiers 161/361...
Training classifiers 162/361...
Training classifiers 163/361...
Training classifiers 164/361...
Training classifiers 165/361...
Training classifiers 166/361...
Training classifiers 167/361...
Training classifiers 168/361...
Training classifiers 169/361...
Training classifiers 170/361...
Training classifiers 171/361...
Training classifiers 172/361...
Training classifiers 173/361...
Training classifiers 174/361...
Training classifiers 175/361...
Training classifiers 176/361...
Training classifiers 177/361...
Training classifiers 178/361...
Training classifiers 179/361...
Training classifiers 180/361...
Training classifiers 181/361...
Training classifiers 182/361...
Training classifiers 183/361...
Training classifiers 184/361...
Training classifiers 185/361...
Training classifiers 186/361...
Training classifiers 187/361...
Training classifiers 188/361...
Training classifiers 189/361...
Training classifiers 190/361...
Training classifiers 191/361...
Training classifiers 192/361...
Training classifiers 193/361...
Training classifiers 194/361...
Training classifiers 195/361...
Training classifiers 196/361...
Training classifiers 197/361...
Training classifiers 198/361...
Training classifiers 199/361...
Training classifiers 200/361...
Training classifiers 201/361...
Training classifiers 202/361...
Training classifiers 203/361...
Training classifiers 204/361...
Training classifiers 205/361...
Training classifiers 206/361...
Training classifiers 207/361...
Training classifiers 208/361...
Training classifiers 209/361...
Training classifiers 210/361...
Training classifiers 211/361...
Training classifiers 212/361...
Training classifiers 213/361...
Training classifiers 214/361...
Training classifiers 215/361...
Training classifiers 216/361...
Training classifiers 217/361...
Training classifiers 218/361...
Training classifiers 219/361...
Training classifiers 220/361...
Training classifiers 221/361...
Training classifiers 222/361...
Training classifiers 223/361...
Training classifiers 224/361...
Training classifiers 225/361...
Training classifiers 226/361...
Training classifiers 227/361...
Training classifiers 228/361...
Training classifiers 229/361...
Training classifiers 230/361...
Training classifiers 231/361...
Training classifiers 232/361...
Training classifiers 233/361...
Training classifiers 234/361...
Training classifiers 235/361...
Training classifiers 236/361...
Training classifiers 237/361...
Training classifiers 238/361...
Training classifiers 239/361...
Training classifiers 240/361...
Training classifiers 241/361...
Training classifiers 242/361...
Training classifiers 243/361...
Training classifiers 244/361...
Training classifiers 245/361...
Training classifiers 246/361...
Training classifiers 247/361...
Training classifiers 248/361...
Training classifiers 249/361...
Training classifiers 250/361...
Training classifiers 251/361...
Training classifiers 252/361...
Training classifiers 253/361...
Training classifiers 254/361...
Training classifiers 255/361...
Training classifiers 256/361...
Training classifiers 257/361...
Training classifiers 258/361...
Training classifiers 259/361...
Training classifiers 260/361...
Training classifiers 261/361...
Training classifiers 262/361...
Training classifiers 263/361...
Training classifiers 264/361...
Training classifiers 265/361...
Training classifiers 266/361...
Training classifiers 267/361...
Training classifiers 268/361...
Training classifiers 269/361...
Training classifiers 270/361...
Training classifiers 271/361...
Training classifiers 272/361...
Training classifiers 273/361...
Training classifiers 274/361...
Training classifiers 275/361...
Training classifiers 276/361...
Training classifiers 277/361...
Training classifiers 278/361...
Training classifiers 279/361...
Training classifiers 280/361...
Training classifiers 281/361...
Training classifiers 282/361...
Training classifiers 283/361...
Training classifiers 284/361...
Training classifiers 285/361...
Training classifiers 286/361...
Training classifiers 287/361...
Training classifiers 288/361...
Training classifiers 289/361...
Training classifiers 290/361...
Training classifiers 291/361...
Training classifiers 292/361...
Training classifiers 293/361...
Training classifiers 294/361...
Training classifiers 295/361...
Training classifiers 296/361...
Training classifiers 297/361...
Training classifiers 298/361...
Training classifiers 299/361...
Training classifiers 300/361...
Training classifiers 301/361...
Training classifiers 302/361...
Training classifiers 303/361...
Training classifiers 304/361...
Training classifiers 305/361...
Training classifiers 306/361...
Training classifiers 307/361...
Training classifiers 308/361...
Training classifiers 309/361...
Training classifiers 310/361...
Training classifiers 311/361...
Training classifiers 312/361...
Training classifiers 313/361...
Training classifiers 314/361...
Training classifiers 315/361...
Training classifiers 316/361...
Training classifiers 317/361...
Training classifiers 318/361...
Training classifiers 319/361...
Training classifiers 320/361...
Training classifiers 321/361...
Training classifiers 322/361...
Training classifiers 323/361...
Training classifiers 324/361...
Training classifiers 325/361...
Training classifiers 326/361...
Training classifiers 327/361...
Training classifiers 328/361...
Training classifiers 329/361...
Training classifiers 330/361...
Training classifiers 331/361...
Training classifiers 332/361...
Training classifiers 333/361...
Training classifiers 334/361...
Training classifiers 335/361...
Training classifiers 336/361...
Training classifiers 337/361...
Training classifiers 338/361...
Training classifiers 339/361...
Training classifiers 340/361...
Training classifiers 341/361...
Training classifiers 342/361...
Training classifiers 343/361...
Training classifiers 344/361...
Training classifiers 345/361...
Training classifiers 346/361...
Training classifiers 347/361...
Training classifiers 348/361...
Training classifiers 349/361...
Training classifiers 350/361...
Training classifiers 351/361...
Training classifiers 352/361...
Training classifiers 353/361...
Training classifiers 354/361...
Training classifiers 355/361...
Training classifiers 356/361...
Training classifiers 357/361...
Training classifiers 358/361...
Training classifiers 359/361...
Training classifiers 360/361...
Training classifiers 361/361...
Training classifiers 361/361... Done (23.40s).

Here we run the prediction: the decoding itself¶

sys.stderr.write("Calculating scores and outputs...")
t0 = time.time()

y_pred = [clf.predict(X_test) for clf in clfs]
y_pred = np.asarray(y_pred).T


# We need to the multi scale reconstruction
def split_multi_scale(y, y_shape):
    """Split data into 4 original multi_scale images"""
    yw, yh = y_shape

    # Index of original image
    split_index = [yw * yh]
    # Index of large image
    split_index.append(split_index[-1] + (yw - 1) * yh)
    # Index of tall image
    split_index.append(split_index[-1] + yw * (yh - 1))
    # Index of big image
    split_index.append(split_index[-1] + (yw - 1) * (yh - 1))

    # We split according to computed indices
    y_preds = np.split(y, split_index, axis=1)

    # y_pred is the original image
    y_pred = y_preds[0]

    # y_pred_tall is the image with 1x2 patch application. We have to make
    # some calculus to get it back in original shape
    height_tf_i = (np.eye(y_cols) + np.eye(y_cols, k=-1))[
        :, : y_cols - 1
    ] * 0.5
    height_tf_i.flat[0] = 1
    height_tf_i.flat[-1] = 1
    y_pred_tall = [
        np.dot(height_tf_i, np.reshape(m, (yw - 1, yh))).flatten()
        for m in y_preds[1]
    ]
    y_pred_tall = np.asarray(y_pred_tall)

    # y_pred_large is the image with 2x1 patch application. We have to make
    # some calculus to get it back in original shape
    width_tf_i = (np.eye(y_cols) + np.eye(y_cols, k=1))[: y_cols - 1] * 0.5
    width_tf_i.flat[0] = 1
    width_tf_i.flat[-1] = 1
    y_pred_large = [
        np.dot(np.reshape(m, (yw, yh - 1)), width_tf_i).flatten()
        for m in y_preds[2]
    ]
    y_pred_large = np.asarray(y_pred_large)

    # y_pred_big is the image with 2x2 patch application. We use previous
    # matrices to get it back in original shape
    y_pred_big = [
        np.dot(np.reshape(m, (yw - 1, yh - 1)), width_tf_i) for m in y_preds[3]
    ]
    y_pred_big = [
        np.dot(height_tf_i, np.reshape(m, (yw - 1, yh))).flatten()
        for m in y_pred_big
    ]
    y_pred_big = np.asarray(y_pred_big)

    return (y_pred, y_pred_tall, y_pred_large, y_pred_big)


y_pred, y_pred_tall, y_pred_large, y_pred_big = split_multi_scale(
    y_pred, y_shape
)

y_pred = (
    0.25 * y_pred
    + 0.25 * y_pred_tall
    + 0.25 * y_pred_large
    + 0.25 * y_pred_big
)

sys.stderr.write(f" Done ({time.time() - t0:.2f}s).\n")

Calculating scores and outputs... Done (1.74s).

Let us quantify our prediction error¶

from sklearn.metrics import (
    accuracy_score,
    f1_score,
    precision_score,
    recall_score,
)

print("Scores")
print("------")
accuracy_to_print = np.mean(
    [accuracy_score(y_test[:, i], y_pred[:, i] > 0.5) for i in range(100)]
)
print(f"  - Accuracy (percent): {accuracy_to_print:f}")

precision_to_print = np.mean(
    [precision_score(y_test[:, i], y_pred[:, i] > 0.5) for i in range(100)]
)
print(f"  - Precision: {precision_to_print:f}")

recall_to_print = np.mean(
    [
        recall_score(y_test[:, i], y_pred[:, i] > 0.5, zero_division=0)
        for i in range(100)
    ]
)
print(f"  - Recall: {recall_to_print:f}")

f1_score_to_print = np.mean(
    [f1_score(y_test[:, i], y_pred[:, i] > 0.5) for i in range(100)]
)
print(f"  - F1-score: {f1_score_to_print:f}")

Scores
------
  - Accuracy (percent): 0.801208
  - Precision: 0.413969
  - Recall: 0.572388
  - F1-score: 0.465973

And finally, we plot six reconstructed images, to compare with ground truth

from pathlib import Path

from matplotlib import pyplot as plt

output_dir = Path.cwd() / "results" / "plot_miyawaki_reconstruction"
output_dir.mkdir(exist_ok=True, parents=True)
print(f"Output will be saved to: {output_dir}")

for i in range(6):
    j = 10 * i
    fig = plt.figure()
    sp1 = plt.subplot(131)
    sp1.axis("off")
    plt.title("Stimulus")
    sp2 = plt.subplot(132)
    sp2.axis("off")
    plt.title("Reconstruction")
    sp3 = plt.subplot(133)
    sp3.axis("off")
    plt.title("Binarized")
    sp1.imshow(
        np.reshape(y_test[j], (10, 10)),
        cmap=plt.cm.gray,
        interpolation="nearest",
    )
    sp2.imshow(
        np.reshape(y_pred[j], (10, 10)),
        cmap=plt.cm.gray,
        interpolation="nearest",
    )
    sp3.imshow(
        np.reshape(y_pred[j] > 0.5, (10, 10)),
        cmap=plt.cm.gray,
        interpolation="nearest",
    )
    plt.savefig(output_dir / f"miyawaki2008_reconstruction_{int(i)}.png")

show()

Output will be saved to: /home/runner/work/nilearn/nilearn/examples/02_decoding/results/plot_miyawaki_reconstruction

References¶

[1]
Yoichi Miyawaki, Hajime Uchida, Okito Yamashita, Masa-aki Sato, Yusuke Morito, Hiroki C. Tanabe, Norihiro Sadato, and Yukiyasu Kamitani. Visual image reconstruction from human brain activity using a combination of multiscale local image decoders. Neuron, 60(5):915–929, 2008. doi:10.1016/j.neuron.2008.11.004.

Total running time of the script: (0 minutes 40.896 seconds)

Estimated memory usage: 1095 MB

Gallery generated by Sphinx-Gallery