OpenScope’s Vismo Dataset

import warnings
warnings.filterwarnings('ignore')

try:
    from databook_utils.dandi_utils import dandi_download_open
except:
    !git clone https://github.com/AllenInstitute/openscope_databook.git
    %cd openscope_databook
    %pip install -e .
    %cd docs/projects

import os

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch
from scipy import interpolate
from scipy.stats import mannwhitneyu
from math import ceil

%matplotlib inline

session_files = pd.read_csv("../../data/vismo_sessions.csv")
session_files

# Function to convert string representation of a set to an actual set
def parse_set_string(set_string):
    # Remove curly braces, split by commas, strip spaces and single quotes
    return set(item.strip("'") for item in set_string.strip("{}").split(", "))

subjects = session_files.groupby('sub_name').agg(
    n_sessions=('session_id', 'nunique'),  # Count unique session IDs
    stim_types=('stim_types', lambda x: set().union(*x.apply(parse_set_string))),  # Union of all sets of stim_types
    genotype=('sub_genotype', 'unique')
).reset_index()

subjects

n_sessions = len(session_files["session_id"].value_counts())
subjects_info = session_files.groupby(["sub_name", "sub_sex"]).size().reset_index().to_dict()
m_count = len([sex for sex in subjects_info["sub_sex"].values() if sex == "M"])
f_count = len([sex for sex in subjects_info["sub_sex"].values() if sex == "F"])

print("Dandiset Overview:")
print(len(session_files), "files")
print(len(subjects_info["sub_name"]), "subjects", m_count, "males", f_count,"females")

dandiset_id = "001196"
dandi_filepath = "sub-762077/sub-762077_ses-multiplane-ophys-762077-2025-03-06-11-49-34-acq-synced_ophys.nwb"
download_loc = "."
dandi_api_key = os.environ["DANDI_API_KEY"]

# This can sometimes take a while depending on the size of the file
io = dandi_download_open(dandiset_id, dandi_filepath, download_loc, dandi_api_key=dandi_api_key)
nwb = io.read()

# view the contents of the nwb interactively
nwb

print("Subject ID",nwb.subject.subject_id)
for plane_name, plane in nwb.lab_meta_data.items():
	print(f"{plane_name} at depth {plane.imaging_depth}")

plane_items = {}

for plane_name, plane in nwb.processing.items():  # iterate over string keys
    try:
        plane_items[plane_name] = (plane['images']['max_projection'], plane['images']['segmentation_mask_image'])
    except (KeyError, TypeError):
        # Skip containers that don't have 'images' or can't be indexed
        continue

n_planes = len(plane_items)
if n_planes == 0:
    print("No planes with images found.")
else:
    # Set number of columns of pairs
    n_cols = 3  # each column will be a pair of images
    n_rows = ceil(n_planes / n_cols)

    fig, axes = plt.subplots(n_rows, n_cols * 2, figsize=(6 * n_cols, 2.5 * n_rows))
    
    # Ensure axes is always 2D
    if n_rows == 1:
        axes = axes.reshape(1, n_cols * 2)

    # Flatten axes for easy indexing
    axes = axes.flatten()

    for idx, (plane_name, (max_proj, seg_mask)) in enumerate(plane_items.items()):
        # Left image of pair
        ax_left = axes[idx * 2]
        ax_left.imshow(max_proj, cmap='gray')
        ax_left.axis('off')
        ax_left.set_title(f"{plane_name} max projection", fontsize=10)

        # Right image of pair
        ax_right = axes[idx * 2 + 1]
        ax_right.imshow(seg_mask, cmap='gray')
        ax_right.axis('off')
        ax_right.set_title(f"{plane_name} segmentation mask", fontsize=10)

    # Turn off any unused axes
    for ax in axes[n_planes * 2:]:
        ax.axis('off')

    plt.tight_layout()
    plt.show()

selected_plane = "VISp_0"

dff = nwb.processing[selected_plane]["dff_timeseries"]["dff_timeseries"]
dff_trace = np.array(dff.data)
dff_timestamps = np.array(dff.timestamps)

print(dff_trace.shape)
print(dff_timestamps.shape)

avg_dff_trace = np.average(dff_trace, axis=1)

%matplotlib inline
n_rois = dff_trace.shape[1]
plt.imshow(dff_trace.transpose(), extent=[dff_timestamps[0], dff_timestamps[-1], 0, n_rois], aspect='auto', vmin=0, vmax=1, interpolation='None')
# plt.yticks(np.arange(n_rois)+0.5, np.arange(n_rois))
plt.ylabel("ROI #")
plt.xlabel("time (s)")
plt.title("DF/F for all ROIs throughout entire session")
cbar = plt.colorbar()
cbar.set_label('DF/F (%)')

for stim_table_name in nwb.intervals.keys():
	print(stim_table_name)

# extract epoch times from stim table where stimulus rows have a different 'block' than following row
# returns list of epochs, where an epoch is of the form (stimulus name, stimulus block, start time, stop time)
def extract_epochs_from_table(stim_name, stim_table, epochs):
    
    # specify a current epoch stop and start time
    epoch_start = stim_table.start_time[0]
    epoch_stop = stim_table.stop_time[0]

    # for each row, try to extend current epoch stop_time
    for i in range(len(stim_table)):
        this_block = stim_table.stim_block[i]
        # if end of table, end the current epoch
        if i+1 >= len(stim_table):
            epochs.append((stim_name, this_block, epoch_start, epoch_stop))
            break
            
        next_block = stim_table.stim_block[i+1]
        # if next row is the same stim block, push back epoch_stop time
        if next_block == this_block:
            epoch_stop = stim_table.stop_time[i+1]
        # otherwise, end the current epoch, start new epoch
        else:
            epochs.append((stim_name, this_block, epoch_start, epoch_stop))
            epoch_start = stim_table.start_time[i+1]
            epoch_stop = stim_table.stop_time[i+1]
    
    return epochs

def extract_all_epochs(nwb):
    # extract epochs from all valid stimulus tables
    epochs = []
    for stim_name in nwb.intervals.keys():
        print(stim_name)
        stim_table = nwb.intervals[stim_name]
        try:
            epochs = extract_epochs_from_table(stim_name, stim_table, epochs)
        except:
            continue
    return epochs

%matplotlib inline
### make plot of chosen trace over time with colored epoch sections

def plot_trace_over_epochs(trace_arrs, timestamps_arrs, epochs, time_start=None, time_end=None, title=None, trace_labels=None, yaxlabels=None, xlabel=None):
    assert len(trace_arrs) == len(timestamps_arrs), "there must be an equal number of traces and timestamps arrays"
    if trace_labels is not None:
        assert len(trace_arrs) == len(trace_labels), "there must be an equal number of traces and trace labels arrays"
    if yaxlabels is not None:
        assert len(trace_arrs) == len(yaxlabels),  "there must be an equal number of traces and y-axis labels arrays"

    fig, ax = plt.subplots(figsize=(15,5), dpi=300)
    if time_start is None:
        time_start = np.min(np.concatenate(timestamps_arrs))
    if time_end is None:
        time_end = np.max(np.concatenate(timestamps_arrs))

    # filter epochs which aren't at least partially in the time window
    bounded_epochs = {epoch for epoch in epochs if epoch[2] < time_end and epoch[3] > time_start}

    # assign unique color to each stimulus name
    stim_names = list({epoch[0] for epoch in bounded_epochs})
    colors = plt.cm.rainbow(np.linspace(0,1,len(stim_names)))
    stim_color_map = {stim_names[i]:colors[i] for i in range(len(stim_names))}

    key = {}
    
    y_hi = np.max(np.concatenate(trace_arrs)) # change these to manually set height of the plot
    y_lo = np.min(np.concatenate(trace_arrs))
    # draw colored rectangles for each epoch
    for epoch in bounded_epochs:
        stim_name, stim_block, epoch_start, epoch_end = epoch
        color = stim_color_map[stim_name]
        rec = ax.add_patch(mpl.patches.Rectangle((epoch_start, y_lo), epoch_end-epoch_start, 100, alpha=0.3, facecolor=color, antialiased=False))
        key[(stim_name)] = rec

    ax.set_xlim(time_start, time_end)
    ax.set_ylim(y_lo, y_hi)
    if xlabel is not None:
        ax.set_xlabel(xlabel)
    if title is not None:
        ax.set_title(title)
    if yaxlabels is not None:
        ax.set_ylabel("\n".join(yaxlabels))

    for i in range(len(trace_arrs)):
        # next_color = plt.rcParams["axes.prop_cycle"].by_key()["color"][i]
        # this_ax = ax.twinx()
        line = ax.plot(timestamps_arrs[i], trace_arrs[i])[0]
        if trace_labels is not None:
            key[(trace_labels[i])] = line

    print(key)
    fig.legend(key.values(), key.keys(), loc="lower right", bbox_to_anchor=(1.25, 0.25))
    
    return ax

epochs = extract_all_epochs(nwb)

# epochs take the form (stimulus name, stimulus block, start time, stop time)
print("Num epochs:",len(epochs))
epochs.sort(key=lambda x: x[2])
for epoch in epochs:
    print(epoch)

# can set these manually to get a closer look at the timeline
time_start = min(epochs, key=lambda epoch: epoch[1])[1]
time_end = max(epochs, key=lambda epoch:epoch[2])[2]
# time_start = 3000
# time_end = 3100

# can set this to change what trace is displayed alongside epochs
display_trace = avg_dff_trace * 100 # to yield percentage
# unit_idx = 30
# display_trace = dff_trace[:,unit_idx] * 100

# adjust these to zoom in to a narrow slice of time
start, stop = None, None
ax = plot_trace_over_epochs([display_trace], [dff_timestamps], epochs, start, stop, "Average DFF and running during a period of sequential session", ['avg dff'], ["Average DF/F"], "Time (s)")
plt.tight_layout()
plt.show()

def interpolate_dff(dff_trace, timestamps, interp_hz):
    """Interpolate dF/F traces onto a regular time grid at interp_hz."""
    time_axis = np.arange(timestamps[0], timestamps[-1], step=1.0 / interp_hz)
    interp_dff = []

    # interpolate channel by channel to save RAM
    for channel in range(dff_trace.shape[1]):
        f = interpolate.interp1d(timestamps, dff_trace[:, channel], axis=0, kind="nearest", fill_value="extrapolate")
        interp_dff.append(f(time_axis))
    interp_dff = np.asarray(interp_dff).T

    return interp_dff, time_axis


def get_dff(nwb, plane, interpolating=True):
    """Return dF/F and timestamps for one imaging plane."""
    series = nwb.processing[plane].data_interfaces["dff_timeseries"].roi_response_series[
        "dff_timeseries"
    ]
    timestamps = series.timestamps[:]
    dff = series.data[:]

    if interpolating:
        interp_hz = nwb.imaging_planes[plane].imaging_rate
        dff, timestamps = interpolate_dff(dff, timestamps, interp_hz)

    return dff, timestamps

plane = "VISp_0"
interpolating = True

dff, timestamps = get_dff(nwb, plane, interpolating=interpolating)
print(dff.shape)

# obtain the stimulus presentation times for the random dot kinematograms
rdkCir_presentations = nwb.intervals["rdkCircle_presentations"].to_dataframe()
rdkSqr_presentations = nwb.intervals["rdkSqr_presentations"].to_dataframe()
rdk_presentations = pd.merge(rdkCir_presentations, rdkSqr_presentations, on=['start_time', 'stop_time'], how='outer', suffixes=('_Circle', '_Sqr'))

full_field_rdk_presentations = rdk_presentations[
    (rdk_presentations['opacity_Circle'] == '0.0') & 
    (rdk_presentations['opacity_Sqr'] == '1.0') & 
    (rdk_presentations['FieldCoherence_Sqr'] == '1.0')
]

full_field_rdk_presentations.head()

full_field_rdk_presentations["Dir_Sqr"].unique()

onset_delay = 0.2
window_pre = 0.5

resp_windows_by_dir = {}
for dir in np.sort(full_field_rdk_presentations['Dir_Sqr'].astype(float).unique()):
    dir_df = full_field_rdk_presentations[full_field_rdk_presentations['Dir_Sqr'] == str(dir)]

    resp_windows = []
    for idx, row in dir_df.iterrows():
        window_start_idx, window_end_idx = np.searchsorted(timestamps, row['start_time'] + onset_delay - window_pre), np.searchsorted(timestamps, row['stop_time'] + onset_delay)
        resp_window = dff[window_start_idx:window_end_idx, :].transpose(1, 0) # window size x n_rois
        resp_windows.append(resp_window)
    
    min_len = min([rw.shape[1] for rw in resp_windows]) # take the minimum length of all response windows to avoid shape mismatch
    resp_windows = np.stack([rw[:, :min_len] for rw in resp_windows], axis=0)
    resp_windows_by_dir[dir] = resp_windows

def show_dff_response(ax, dff, window_start_time, window_end_time, aspect="auto", vmin=None, vmax=None, yticklabels=[], skipticks=1, xlabel="Time (s)", ylabel="ROI", cbar=True, cbar_label=None):
    if len(dff) == 0:
        print("Input data has length 0; Nothing to display")
        return

    img = ax.imshow(dff, aspect=aspect, extent=[window_start_time, window_end_time, 0, len(dff)], vmin=vmin, vmax=vmax)
    if cbar:
        ax.colorbar(img, shrink=0.5, label=cbar_label)

    ax.plot([0,0],[0, len(dff)], ":", color="white", linewidth=1.0)

    if len(yticklabels) != 0:
        ax.set_yticks(range(len(yticklabels)))
        ax.set_yticklabels(yticklabels, fontsize=8)

        n_ticks = len(yticklabels[::skipticks])
        ax.yaxis.set_major_locator(plt.MaxNLocator(n_ticks))

    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)

def show_rois_by_dir(resp_windows_by_dir, window_start_time, window_end_time, vmin=None, vmax=None, cbar_label=None, roi_idxs=None, suptitle=None
):
    """Plot rois x conditions with one colorbar per row."""
    dirs = list(resp_windows_by_dir.keys())
    n_rois = resp_windows_by_dir[dirs[0]].shape[1]
    n_dirs = len(dirs)

    if roi_idxs is None:
        roi_idxs = list(range(n_rois))
    n_rois_to_plot = len(roi_idxs)

    fig, axes = plt.subplots(n_rois_to_plot, n_dirs, figsize=(2*n_dirs, 1.5*n_rois_to_plot), squeeze=False, layout="constrained")
    low_q, high_q = 2, 98  # or 5, 95

    if vmin is None or vmax is None:
        row_vmin, row_vmax = [], []
        for r_idx_data in roi_idxs:
            vals = np.concatenate([resp_windows_by_dir[dir][:, r_idx_data, :].ravel() for dir in dirs])
            row_vmin.append(np.nanpercentile(vals, low_q))
            row_vmax.append(np.nanpercentile(vals, high_q))

        if row_vmax[-1] <= row_vmin[-1]:
            row_vmax[-1] = row_vmin[-1] + 1e-6    
    else:    
        row_vmin = [vmin] * n_rois_to_plot
        row_vmax = [vmax] * n_rois_to_plot

    for c_idx, dir in enumerate(dirs):
        for r_idx_plot, r_idx_data in enumerate(roi_idxs):
            ax = axes[r_idx_plot, c_idx]
            dff = resp_windows_by_dir[dir][:, r_idx_data, :]

            this_vmin = row_vmin[r_idx_plot]
            this_vmax = row_vmax[r_idx_plot]

            show_dff_response(
                ax, dff,
                window_start_time, window_end_time,
                xlabel="Time (frames)" if r_idx_plot == n_rois_to_plot-1 else None,
                ylabel="Trials" if c_idx == 0 else None,
                cbar=False,
                vmin=this_vmin, vmax=this_vmax,
            )

            if r_idx_plot == 0:
                ax.set_title(f"{float(dir):.0f}°", fontsize=11, pad=10)
            if c_idx == 0:
                ax.set_ylabel(f"ROI {r_idx_data}\nTrials")

    for r_idx_plot in range(n_rois_to_plot):
        norm = mpl.colors.Normalize(vmin=row_vmin[r_idx_plot], vmax=row_vmax[r_idx_plot])
        fig.colorbar(
            mpl.cm.ScalarMappable(norm=norm),
            ax=axes[r_idx_plot, :],
            orientation='vertical',
            shrink=0.6, pad=0.02,
            label=cbar_label,
        )

    if suptitle is not None:
        fig.suptitle(suptitle, fontsize=12)

    return fig, axes

# simple heuristic to select cells to plot
def select_top_rois_tuning_ampl(windows, pre_n=5, on_n=10, top_k=10):
    """Score ROIs by direction-tuning prominence and return the top_k indices.
    Uses (on - pre) responses per direction; score = peak direction minus mean of other directions.
    """
    per_dir_means = []  # list of (cells,) arrays

    for w in windows:
        n_time = w.shape[-1]
        pre = w[:, :, :min(pre_n, n_time // 2)].mean(-1)       
        on  = w[:, :, -min(on_n, n_time // 2):].mean(-1)      
        delta_mean = (on - pre).mean(axis=0)                   
        per_dir_means.append(delta_mean)

    mean_dir = np.stack(per_dir_means, axis=0)                 

    peak = mean_dir.max(axis=0)                                
    peak_dir = mean_dir.argmax(axis=0)                         

    mean_other = (mean_dir.sum(axis=0) - peak) / (mean_dir.shape[0] - 1)
    prominence = peak - mean_other                             
    # require positive peak response
    prominence[peak <= 0] = -np.inf

    idx = np.argsort(prominence)[-top_k:][::-1]
    return idx, prominence

roi_idxs, _ = select_top_rois_tuning_ampl(resp_windows_by_dir.values(), top_k=10)
print("Selected ROIs:", roi_idxs)

fig, axes = show_rois_by_dir(
    resp_windows_by_dir=resp_windows_by_dir,
    window_start_time=-5,             
    window_end_time=10,                
    vmin=None,                          
    vmax=None,
    cbar_label="dF/F",
    roi_idxs=roi_idxs,
    suptitle="Responses by direction",
)

def plot_tuning_curves(
    resp_windows_by_dir,
    roi_indices,
    on_n_frames=5,
    n_rows=2,
    n_cols=3,
):
    """Plot tuning curves with trial scatter and per-ROI 98th percentile clipping."""
    directions_deg = list(resp_windows_by_dir.keys())
    roi_indices = list(roi_indices)

    n_cells = resp_windows_by_dir[directions_deg[0]].shape[1]
    n_plots = min(len(roi_indices), n_rows * n_cols)

    fig, axes = plt.subplots(n_rows, n_cols, figsize=(3 * n_cols, 3 * n_rows), sharex=True, squeeze=False)
    axes = axes.ravel()
    for ax in axes[n_plots:]:
        ax.set_visible(False)

    color_mean = "black"
    color_dots = "gray"

    for plot_idx, roi_idx in enumerate(roi_indices[:n_plots]):
        ax = axes[plot_idx]
        if roi_idx >= n_cells:
            ax.set_visible(False)
            continue

        means_per_dir = []
        all_vals = []
        dir_maxes = []
        dir_vals = []

        for dir in directions_deg:
            w = resp_windows_by_dir[dir]
            if roi_idx >= w.shape[1]:
                continue
            n_time = w.shape[-1]
            on_n = min(on_n_frames, n_time)
            trial_vals = w[:, roi_idx, -on_n:].mean(axis=-1)
            means_per_dir.append(trial_vals.mean())
            all_vals.append(trial_vals)
            dir_maxes.append(trial_vals.max())
            dir_vals.append(dir)
            ax.scatter(np.full_like(trial_vals, dir_vals[-1], dtype=float), trial_vals, color=color_dots, s=15, alpha=0.3)

        if all_vals:
            all_vals = np.concatenate(all_vals)
            y_max = np.nanpercentile(all_vals, 98)
            ax.set_ylim(top=y_max * 1.05)
            for dv, dmax in zip(dir_vals, dir_maxes):
                if dmax > y_max:
                    ax.scatter(dv, y_max, marker="^", color=color_dots, s=25)

        means_per_dir = np.asarray(means_per_dir)
        ax.plot(directions_deg, means_per_dir, "o-", color=color_mean, linewidth=2)
        ax.axhline(0, color="gray", linestyle=":", linewidth=0.8)
        ax.set_title(f"ROI {roi_idx}")
        ax.set_xticks(list(resp_windows_by_dir.keys()))
        if plot_idx // n_cols == n_rows - 1:
            ax.set_xlabel("Direction (°)")
        if plot_idx % n_cols == 0:
            ax.set_ylabel("Response (dF/F)")

    plt.tight_layout()
    return fig

n_cols = 5
n_rows = roi_idxs.size // n_cols
fig = plot_tuning_curves(resp_windows_by_dir, roi_idxs, on_n_frames=10, n_rows=n_rows, n_cols=n_cols)

# compute the pre-stimulus and during-stimulus responses for each stimulus presentation
pre_stim_start_time = -5 # in frames
indicator_delay = 2 # in frames

def compute_pre_and_during_stimulus_response(row):
    assert pre_stim_start_time < 0, "pre_stim_start_time must be negative"
    stim_start_idx, stim_stop_idx = np.searchsorted(dff_timestamps, row['start_time']) + indicator_delay, np.searchsorted(dff_timestamps, row['stop_time']) + indicator_delay
    pre_stim_resp = dff_trace[stim_start_idx+pre_stim_start_time:stim_start_idx].mean(axis=0)
    stim_resp = dff_trace[stim_start_idx:stim_stop_idx].mean(axis=0)
    return pd.Series([pre_stim_resp, stim_resp])

rdk_presentations[['pre_stim_resp', 'stim_resp']] = rdk_presentations.apply(compute_pre_and_during_stimulus_response, axis=1)

# use the wilcoxon rank-sum test to compare the pre-stimulus responses and during-stimulus responses
pre_stim_resp = np.vstack(rdk_presentations['pre_stim_resp'].values)
stim_resp = np.vstack(rdk_presentations['stim_resp'].values)

p_vals = np.array([
    mannwhitneyu(stim_resp[:, i], pre_stim_resp[:, i], alternative='greater').pvalue
    for i in range(stim_resp.shape[1])
])

# Define responsiveness
responsive = p_vals < 0.05
cmap = ListedColormap(['tab:red', 'tab:green'])

fig, ax = plt.subplots(figsize=(10, 1.5))

ax.imshow(responsive[np.newaxis, :], aspect='auto', cmap=cmap)
ax.set_yticks([])
ax.set_xlabel('ROI index')
ax.set_title('ROI Responsiveness (Wilcoxon rank-sum test)')
for spine in ax.spines.values():
    spine.set_visible(False)
legend_elements = [
    Patch(facecolor='tab:green', edgecolor='none', label='Responsive (p < 0.05)'),
    Patch(facecolor='tab:red', edgecolor='none', label='Not responsive (p ≥ 0.05)')
]
ax.legend(handles=legend_elements, loc='upper right', frameon=False, bbox_to_anchor=(1.25, 0.05))

plt.tight_layout()
plt.show()

rf_block_presentations = nwb.intervals["receptive_field_block_presentations"].to_dataframe()

def compute_stim_response(row):
    stim_start_idx = np.searchsorted(dff_timestamps, row['start_time'] + onset_delay)
    stim_stop_idx = np.searchsorted(dff_timestamps, row['stop_time'] + onset_delay)
    
    return dff_trace[stim_start_idx:stim_stop_idx].mean(axis=0)

rf_block_presentations['resp'] = rf_block_presentations.apply(compute_stim_response, axis=1)

resps = np.vstack(rf_block_presentations['resp'].values)
nrois = resps.shape[1]

# create receptive field maps
unique_x = np.sort(np.unique(rf_block_presentations['x_position'].astype(float)))
unique_y = np.sort(np.unique(rf_block_presentations['y_position'].astype(float)))

receptive_fields = np.zeros((nrois, len(unique_y), len(unique_x)))
stim_counts = np.zeros((len(unique_y), len(unique_x)))
for i, row in rf_block_presentations.iterrows():
    x_pos = float(row['x_position'])
    y_pos = float(row['y_position'])

    x_idx = np.argmin(np.abs(unique_x - x_pos))
    y_idx = np.argmin(np.abs(unique_y - y_pos))

    resp = row['resp']
    receptive_fields[:, y_idx, x_idx] += resp  # No need for [:, None]
    stim_counts[y_idx, x_idx] += 1

receptive_fields /= stim_counts[None, :, :]

nrows, ncols = 2, 5
fig, axes = plt.subplots(nrows, ncols, figsize=(ncols*5, nrows*4), constrained_layout=True)
axes = axes.ravel()

for i in range(nrows*ncols):
    ax = axes[i]
    roi_idx = i
    im = ax.imshow(
        receptive_fields[roi_idx],
        extent=[float(unique_x.min()), float(unique_x.max()), float(unique_y.min()), float(unique_y.max())],
        aspect='auto',
        origin='upper'
    )
    ax.set_title(f"ROI {roi_idx}", fontsize=15)
    ax.set_xlabel("X position (deg)", fontsize=12)
    ax.set_ylabel("Y position (deg)", fontsize=12)
    cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
    cbar.set_label('Response', fontsize=12)


plt.show()