Decoding and Identifying Natural Images from fMRI

Ben Oliver. Final project for CSCI 591 at the University of Montana (see kay_pitch.ipynb for the original pitch)

Data: Kay, Naselaris, Prenger & Gallant (2008), Identifying natural images from human brain activity, Nature 452.

Summary

We ask how much of a natural image can be recovered from fMRI responses in visual cortex, and which visual areas carry that information. Two approaches are compared on the dataset's 120 held-out test images:

1. Decoding (brain to pixels). Ridge regression from voxel responses to the image's principal components. Reconstructions recover only the coarsest layout of a scene, mostly whether its top or bottom is brighter. Using early visual cortex (V1 to V3), the decoder can tell the true image from another test image 66% of the time (chance 50%), and picks the exact image out of 120 only 2.5% of the time (chance 0.8%).
2. Encoding (pixels to brain), then identification. Following Kay et al., we model each voxel's response as a weighted sum of Gabor energy features, then identify which image was seen by asking which candidate image's predicted brain response best matches the measured one. Using the 500 best-predicted voxels this identifies the correct image out of 120 78% of the time (93 of 120), and wins the pairwise comparison 99% of the time.

In both approaches, V1 and V2 carry the most image information and accuracy drops sharply in V3 and beyond (V3A, V3B, V4 and lateral occipital cortex). That is what we expect from areas whose voxels are tuned to local orientation and spatial frequency, and it shows why simple linear models do far better at identifying an image than at reconstructing it.

1. Data

The dataset has one subject's responses to 1,750 training images and 120 test images (each test image was shown 13 times and the responses averaged). Images are 128 x 128 grayscale. Responses are z-scored, voxel-wise BOLD amplitudes for 8,428 voxels, each assigned to a visual area.

import os
import requests
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter
from skimage.metrics import structural_similarity
from sklearn.decomposition import PCA
from sklearn.linear_model import RidgeCV

plt.rcParams.update({"figure.dpi": 110, "axes.spines.top": False, "axes.spines.right": False})
rng = np.random.default_rng(0)

# download the data if it is not already present
for fname, url in [("kay_images.npz", "https://osf.io/ymnjv/download")]:
    if not os.path.isfile(fname):
        r = requests.get(url, timeout=600)
        r.raise_for_status()
        with open(fname, "wb") as f:
            f.write(r.content)

with np.load("kay_images.npz") as d:
    dat = dict(**d)

roi_names = [str(r) for r in dat["roi_names"]]
X_tr = dat["responses"].astype(np.float32)        # (1750, 8428) training responses
X_te = dat["responses_test"].astype(np.float32)   # (120, 8428) test responses
S_tr, S_te = dat["stimuli"], dat["stimuli_test"]  # images
Y_tr, Y_te = S_tr.reshape(len(S_tr), -1), S_te.reshape(len(S_te), -1)
n_test = len(S_te)

print("training:", X_tr.shape, "test:", X_te.shape)
print({n: int((dat["roi"] == i).sum()) for i, n in enumerate(roi_names) if i})

training: (1750, 8428) test: (120, 8428)
{'V1': 1294, 'V2': 2083, 'V3': 1790, 'V3A': 484, 'V3B': 314, 'V4': 1535, 'LatOcc': 928}

fig, axs = plt.subplots(1, 6, figsize=(12, 2.4))
for ax, i in zip(axs, rng.choice(n_test, 6, replace=False)):
    ax.imshow(S_te[i], cmap="gray")
    ax.axis("off")
fig.suptitle("Example test images")
plt.show()

2. How we measure success

• Pixel correlation. Every image in this dataset sits inside the same circular mask, so even the average training image correlates with each test image (r about 0.2). To measure what the brain data adds, we subtract the mean training image from both the reconstruction and the truth before correlating. We call this the residual correlation.
• SSIM (structural similarity) was in the original plan. We report it once below, but the mean image alone scores the same as every decoder, so it does not separate models here.
• Pairwise identification. For each test image, is its reconstruction (or predicted response) closer to the true image than to each of the other 119? Chance is 50%.
• Top-1 identification. Out of all 120 test images, is the best match the right one? Chance is 1/120 = 0.8%.

def zrows(a):
    a = a - a.mean(1, keepdims=True)
    return a / np.linalg.norm(a, axis=1, keepdims=True)

def identification(A, B):
    # C[i, j] = correlation between A_i and B_j; row i is correct if C[i, i] is the largest
    C = zrows(A) @ zrows(B).T
    top1 = np.mean(C.argmax(1) == np.arange(len(C)))
    off = ~np.eye(len(C), dtype=bool)
    pairwise = np.mean((np.diag(C)[:, None] > C)[off])
    return top1, pairwise, np.diag(C)

def mean_ssim(P, T):
    return np.mean([structural_similarity(p.reshape(128, 128), t.reshape(128, 128), data_range=t.max() - t.min())
                    for p, t in zip(P, T)])

mean_img = Y_tr.mean(0)
baseline = np.repeat(mean_img[None], n_test, 0)
print("mean training image as the 'reconstruction': r = %.3f, SSIM = %.3f"
      % (identification(baseline, Y_te)[2].mean(), mean_ssim(baseline, Y_te)))

mean training image as the 'reconstruction': r = 0.215, SSIM = 0.210

3. Decoding: from brain activity to pixels

We compress the training images with PCA (128 components) and fit a ridge regression from voxel responses to each component's score, choosing a separate regularization strength for each component by efficient leave-one-out cross-validation on the training set. The predicted scores are mapped back to pixels. Predicting PCA scores instead of raw pixels gives slightly better reconstructions and is much faster. We fit one decoder per visual area, one for early visual cortex combined (V1 to V3), and one on all voxels.

pca = PCA(128, random_state=0).fit(Y_tr)
Z_tr = pca.transform(Y_tr)
alphas = np.logspace(1, 6, 11)

rois = {n: dat["roi"] == i for i, n in enumerate(roi_names) if i}
rois["V1-V3"] = np.isin(dat["roi"], [1, 2, 3])
rois["All"] = np.ones(X_tr.shape[1], bool)

decode = {}
for name, m in rois.items():
    model = RidgeCV(alphas=alphas, alpha_per_target=True).fit(X_tr[:, m], Z_tr)
    P = pca.inverse_transform(model.predict(X_te[:, m]))
    top1, pair, r = identification(P - mean_img, Y_te - mean_img)
    decode[name] = dict(P=P, r=r, top1=top1, pair=pair)
    print(f"{name:7s} voxels={m.sum():5d}  residual r={r.mean():.3f}  pairwise={pair:.1%}  top-1={top1:.1%}")

print("\nSSIM, V1-V3 decoder: %.3f (mean image: %.3f)" % (mean_ssim(decode["V1-V3"]["P"], Y_te), mean_ssim(baseline, Y_te)))

V1      voxels= 1294  residual r=0.151  pairwise=63.8%  top-1=3.3%

V2      voxels= 2083  residual r=0.148  pairwise=62.3%  top-1=2.5%

V3      voxels= 1790  residual r=0.131  pairwise=58.2%  top-1=1.7%
V3A     voxels=  484  residual r=0.053  pairwise=55.6%  top-1=2.5%

V3B     voxels=  314  residual r=0.120  pairwise=59.6%  top-1=0.8%

V4      voxels= 1535  residual r=0.130  pairwise=58.6%  top-1=2.5%

LatOcc  voxels=  928  residual r=0.077  pairwise=55.0%  top-1=0.0%

V1-V3   voxels= 5167  residual r=0.176  pairwise=65.7%  top-1=2.5%

All     voxels= 8428  residual r=0.198  pairwise=66.8%  top-1=3.3%

SSIM, V1-V3 decoder: 0.214 (mean image: 0.210)

best = decode["V1-V3"]
order = np.argsort(best["r"])[::-1]
show = list(order[:4]) + list(order[n_test // 2 - 1:n_test // 2 + 1]) + list(order[-2:])

fig, axs = plt.subplots(2, len(show), figsize=(14, 3.8))
for k, i in enumerate(show):
    axs[0, k].imshow(S_te[i], cmap="gray")
    axs[1, k].imshow(best["P"][i].reshape(128, 128), cmap="gray")
    axs[1, k].set_title(f"r = {best['r'][i]:.2f}", fontsize=9)
    for ax in axs[:, k]:
        ax.axis("off")
axs[0, 0].text(-12, 64, "true", rotation=90, va="center", ha="right")
axs[1, 0].text(-12, 64, "decoded", rotation=90, va="center", ha="right")
fig.suptitle("Decoded from V1-V3: four best, two middle, two worst (residual correlation)")
plt.show()

The reconstructions all look alike: a blurry disk that is brighter at the top or at the bottom. The best-scoring images are scenes with a bright sky over darker ground, and the worst are scenes lit the other way. The decoder recovers that vertical gradient but none of the edges or texture. One way to check this: compare the reconstructions with blurred versions of the true images. If the decoder only captures coarse structure, it should match the blurred images better than the sharp ones.

sigmas = [0, 1, 2, 4, 6, 8, 12]
blur_r = []
for s in sigmas:
    T = np.stack([gaussian_filter(t.reshape(128, 128), s).ravel() for t in Y_te]) if s else Y_te
    blur_r.append(identification(best["P"] - mean_img, T - mean_img)[2].mean())

fig, ax = plt.subplots(figsize=(5, 3.2))
ax.plot(sigmas, blur_r, marker="o")
ax.set_xlabel("blur applied to true image (Gaussian sigma, pixels)")
ax.set_ylabel("residual correlation")
ax.set_title("Decoded images match blurred truth better")
plt.show()

Correlation keeps rising as the true image is blurred, so what the decoder recovers is mostly low spatial frequency. That is the "blurry reconstruction" problem we expected in the pitch: a linear map from about 5,000 noisy voxels, trained on 1,750 images, averages away fine detail.

4. Encoding, then identification (the Kay et al. approach)

Instead of predicting the image from the brain, we predict the brain from the image, where the model can use what we know about early vision:

1. Features. Each image is filtered with a bank of Gabor filters (5 spatial frequencies from 2 to 32 cycles per image, 8 orientations). We keep the phase-invariant energy of each filter and average it over a grid of positions that gets finer for higher frequencies, giving 6,784 features per image. This is a simplified version of the Gabor wavelet pyramid in Kay et al.
2. Voxel models. A separate ridge regression for each voxel predicts its response from the features, with its own regularization strength.
3. Voxel selection. We hold out 10% of the training images, measure how well each voxel's model predicts them, and keep the best voxels. The test set plays no part in this choice.
4. Identification. For a measured test response, we predict the response to every one of the 120 candidate images and pick the candidate whose prediction correlates best with what was measured.

def gabor_energy(imgs, freqs=(2, 4, 8, 16, 32), n_orient=8, max_grid=16, bw=0.35):
    # A one-sided Gaussian in the Fourier domain is a complex (quadrature-pair) Gabor in space,
    # so the magnitude of the filtered image is the phase-invariant Gabor energy.
    n, h, w = imgs.shape
    F = np.fft.fft2(imgs)
    fy, fx = np.meshgrid(np.fft.fftfreq(h) * h, np.fft.fftfreq(w) * w, indexing="ij")  # cycles per image
    feats = []
    for f in freqs:
        g = int(min(2 * f, max_grid))  # pooling grid: g x g positions
        for th in np.arange(n_orient) * np.pi / n_orient:
            G = np.exp(-((fx - f * np.cos(th)) ** 2 + (fy - f * np.sin(th)) ** 2) / (2 * (bw * f) ** 2))
            E = np.abs(np.fft.ifft2(F * G)).reshape(n, g, h // g, g, w // g).mean((2, 4))
            feats.append(np.sqrt(E.reshape(n, -1)))  # square root compresses large energies
    return np.concatenate(feats, 1).astype(np.float32)

Phi = gabor_energy(np.concatenate([S_tr, S_te]))
Phi_tr, Phi_te = Phi[:len(S_tr)], Phi[len(S_tr):]
mu, sd = Phi_tr.mean(0), Phi_tr.std(0) + 1e-8
Phi_tr, Phi_te = (Phi_tr - mu) / sd, (Phi_te - mu) / sd
print("features per image:", Phi_tr.shape[1])

features per image: 6784

enc_alphas = np.logspace(2, 6, 9)

# voxel selection on a split of the training set only
val = rng.permutation(len(Phi_tr))[:175]
fit = np.setdiff1d(np.arange(len(Phi_tr)), val)
sel_model = RidgeCV(alphas=enc_alphas, alpha_per_target=True).fit(Phi_tr[fit], X_tr[fit])
zcols = lambda a: (a - a.mean(0)) / a.std(0)
vox_r = (zcols(sel_model.predict(Phi_tr[val])) * zcols(X_tr[val])).mean(0)

# final voxel models on all training images; predicted response to each test image
enc = RidgeCV(alphas=enc_alphas, alpha_per_target=True).fit(Phi_tr, X_tr)
P_te = enc.predict(Phi_te)

def identify(vox):
    top1, pair, _ = identification(X_te[:, vox], P_te[:, vox])
    return top1, pair

top500 = np.argsort(vox_r)[::-1][:500]
top1, pair = identify(top500)
print(f"500 best voxels: top-1 = {top1:.1%} (chance {1/n_test:.1%}), pairwise = {pair:.1%}")
print("visual areas of those 500 voxels:", {n: int((dat['roi'][top500] == i).sum()) for i, n in enumerate(roi_names) if i})

500 best voxels: top-1 = 77.5% (chance 0.8%), pairwise = 98.8%
visual areas of those 500 voxels: {'V1': 207, 'V2': 209, 'V3': 48, 'V3A': 0, 'V3B': 4, 'V4': 24, 'LatOcc': 8}

areas = [n for n in roi_names if n != "Other"]
fig, axs = plt.subplots(1, 2, figsize=(12, 3.6))

axs[0].boxplot([vox_r[dat["roi"] == roi_names.index(n)] for n in areas], showfliers=False)
axs[0].set_xticks(range(1, len(areas) + 1), areas)
axs[0].axhline(0, color="gray", lw=0.8)
axs[0].set_ylabel("prediction r on held-out training images")
axs[0].set_title("How well the Gabor model predicts each voxel")

per_area = {}
for n in areas:
    idx = np.where(dat["roi"] == roi_names.index(n))[0]
    per_area[n] = identify(idx[np.argsort(vox_r[idx])[::-1][:250]])
axs[1].bar(areas, [per_area[n][0] for n in areas], color="tab:blue")
axs[1].axhline(1 / n_test, color="gray", ls="--", lw=0.8, label="chance")
axs[1].set_ylabel("top-1 identification (of 120)")
axs[1].set_title("Identification from the 250 best voxels in each area")
axs[1].legend()
plt.tight_layout()
plt.show()

for n in areas:
    print(f"{n:7s} top-1 = {per_area[n][0]:5.1%}   pairwise = {per_area[n][1]:5.1%}")

V1      top-1 = 74.2%   pairwise = 98.7%
V2      top-1 = 62.5%   pairwise = 96.4%
V3      top-1 = 30.0%   pairwise = 90.9%
V3A     top-1 =  9.2%   pairwise = 69.9%
V3B     top-1 =  7.5%   pairwise = 74.1%
V4      top-1 = 15.0%   pairwise = 84.3%
LatOcc  top-1 =  6.7%   pairwise = 77.9%

Ns = [50, 100, 250, 500, 1000, 2000, 4000]
acc = [identify(np.argsort(vox_r)[::-1][:N])[0] for N in Ns]
fig, ax = plt.subplots(figsize=(5, 3.2))
ax.semilogx(Ns, acc, marker="o")
ax.set_xlabel("number of voxels (best first)")
ax.set_ylabel("top-1 identification (of 120)")
ax.set_title("Identification vs. number of voxels")
plt.show()

Accuracy rises quickly over the first few hundred voxels, peaks around 500, and then falls as poorly predicted voxels add more noise than signal. We report 500 voxels as the headline number because that was fixed in advance (it is also the number Kay et al. used), not picked from this curve.

5. Comparison

Approach	Voxels	Pairwise (chance 50%)	Top-1 of 120 (chance 0.8%)
Decoding to pixels, V1-V3	5,167	65.7%	2.5%
Decoding to pixels, all voxels	8,428	66.8%	3.3%
Gabor encoding model + identification	best 500	98.8%	77.5%

Kay et al. (2008) report 92% and 72% top-1 on 120 images for their two subjects with a fuller Gabor pyramid and receptive-field model, so our simplified version (77.5%) lands between them.

6. Discussion

Why identification beats reconstruction. Decoding has to produce all 16,384 pixels from noisy voxels, and ridge regression handles that uncertainty by shrinking toward the average image, which removes detail. Identification only has to rank 120 candidates, and the encoding model can use what we know about early visual cortex: voxels respond to local orientation and spatial frequency energy, and not to the sign of each pixel. That knowledge is built into the Gabor features and is exactly what a pixel-space decoder has to learn from scratch.

Which areas carry the information. Both approaches put V1 and V2 well ahead of everything else. With the encoding model, the 250 best V1 voxels alone identify 74% of images and V2 63%, while V3 manages 30% and the higher areas 7 to 15%. Of the 500 voxels chosen from the whole brain, 416 are in V1 or V2. These higher areas respond to more abstract properties (shape, objects, categories) that Gabor energy and pixel intensity do not capture, so a linear model on these features cannot use them. This matches the literature review: low-level models recover layout and texture from early visual cortex, and recovering semantic content from higher areas needs richer features such as deep network activations.

Limitations. This is one subject. Our Gabor bank is coarser than Kay et al.'s and we did not model each voxel's receptive field explicitly. The decoder is linear. Results could change with a different random split for voxel selection, although the overall pattern does not.

Next steps.

• Decode Gabor or deep network (for example VGG) features instead of pixels, then reconstruct by searching for an image that matches them (Shen et al., 2019).
• Use the encoding model to reconstruct: score a large library of natural images by how well their predicted responses match the measured response, and average the best matches (Naselaris et al., 2009).
• Test whether adding higher areas helps once the features include semantic information, using the WordNet labels that come with the dataset.