Pairwise surface alignment.

In this tutorial, we show how to align surface data from two subjects using a pairwise alignment method. We project the data on the fsaverage5 surface and learn a piecewise mapping from one subject to the other using ward clustering.

We mostly rely on python common packages and on nilearn to handle functional data in a clean fashion.

Retrieve the data

In this example we use the IBC dataset, which include a large number of different contrasts maps for 12 subjects. We download the images for subjects sub-01 and sub-04 (or retrieve them if they were already downloaded) Files is the list of paths for each subjects. df is a dataframe with metadata about each of them.

from fmralign.fetch_example_data import fetch_ibc_subjects_contrasts

files, df, _ = fetch_ibc_subjects_contrasts(["sub-01", "sub-04"])

[get_dataset_dir] Dataset found in /home/runner/nilearn_data/ibc

Project the data on the fsaverage5 surface

We project the data on the fsaverage5 surface, using the fsaverage5 surface template.

from nilearn.datasets import load_fsaverage, load_fsaverage_data
from nilearn.surface import SurfaceImage

fsaverage_meshes = load_fsaverage()


def project_to_surface(img):
    """Util function for loading and projecting volumetric images."""
    surface_image = SurfaceImage.from_volume(
        mesh=fsaverage_meshes["pial"],
        volume_img=img,
    )
    return surface_image


from nilearn.image import concat_imgs

source_train = concat_imgs(
    df[df.subject == "sub-01"][df.acquisition == "ap"].path.values
)
target_train = concat_imgs(
    df[df.subject == "sub-04"][df.acquisition == "ap"].path.values
)

surf_source_train = project_to_surface(source_train)
surf_target_train = project_to_surface(target_train)

/home/runner/work/fmralign/fmralign/examples/plot_surf_alignment.py:53: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  df[df.subject == "sub-01"][df.acquisition == "ap"].path.values
/home/runner/work/fmralign/fmralign/examples/plot_surf_alignment.py:56: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  df[df.subject == "sub-04"][df.acquisition == "ap"].path.values

Fitting the alignment operator

We use the PairwiseAlignment class to learn the alignment operator from one subject to the other. We select the scaled_orthogonal method to compute a rigid piecewise alignment mapping and the ward clustering method to parcellate the cortical surface.

from fmralign.pairwise_alignment import PairwiseAlignment

alignment_estimator = PairwiseAlignment(
    alignment_method="scaled_orthogonal",
    n_pieces=100,
    clustering="ward",
)
# Learn alignment operator from subject 1 to subject 2 on training data
alignment_estimator.fit(surf_source_train, surf_target_train)

/home/runner/work/fmralign/fmralign/fmralign/_utils.py:258: UserWarning: Overriding provided-default estimator parameters with provided masker parameters :
Parameter smoothing_fwhm :
    Masker parameter None - overriding estimator parameter 4.0

  parcellation.fit(images_to_parcel)
/home/runner/work/fmralign/fmralign/.venv/lib/python3.12/site-packages/scipy/sparse/_index.py:210: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil and dok are more efficient.
  self._set_arrayXarray(i, j, x)
/home/runner/work/fmralign/fmralign/.venv/lib/python3.12/site-packages/sklearn/cluster/_agglomerative.py:321: UserWarning: the number of connected components of the connectivity matrix is 2 > 1. Completing it to avoid stopping the tree early.
  connectivity, n_connected_components = _fix_connectivity(
/home/runner/work/fmralign/fmralign/fmralign/_utils.py:189: UserWarning:
 Some parcels are more than 1000 voxels wide it can slow down alignment,especially optimal_transport :
 parcel 17 : 1909 voxels
  warnings.warn(warning)

PairwiseAlignment(alignment_method='scaled_orthogonal',
                  masker=SurfaceMasker(memory=Memory(location=None),
                                       memory_level=0),
                  n_pieces=100)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Plot the computed parcellation

We can retrieve the computed parcellation and plot it on the surface.

from nilearn import plotting

_, clustering_img = alignment_estimator.get_parcellation()


plotting.plot_surf_roi(
    surf_mesh=fsaverage_meshes["pial"],
    roi_map=clustering_img,
    hemi="left",
    view="lateral",
    title="Ward parcellation on the left hemisphere",
)
plotting.show()

/home/runner/work/fmralign/fmralign/examples/plot_surf_alignment.py:91: DeprecationWarning: The `darkness` parameter will be deprecated in release 0.13. We recommend setting `darkness` to None
  plotting.plot_surf_roi(

Projecting the left-out data

Let’s now align a left-out audio contrast from sub-01 to sub-04. We project the data on the surface and apply the learned alignment operator.

surf_audio_source = project_to_surface(
    df[
        (df.subject == "sub-01")
        & (df.condition == "audio_sentence")
        & (df.acquisition == "pa")
    ].path.values
)

surf_audio_target = project_to_surface(
    df[
        (df.subject == "sub-04")
        & (df.condition == "audio_sentence")
        & (df.acquisition == "pa")
    ].path.values
)

surf_aligned = alignment_estimator.transform(surf_audio_source)

Visualizing the alignment in action

We interpolate between the source and aligned images to visualize the alignment process. Notice how the individual idiocyncracies of the source subject are progressively removed.

from copy import deepcopy

import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

fsaverage_sulcal = load_fsaverage_data(
    mesh="fsaverage5",
    data_type="sulcal",
    mesh_type="inflated",
)

plotting_params = {
    "bg_map": fsaverage_sulcal,
    "hemi": "left",
    "view": "lateral",
    "colorbar": True,
    "alpha": 0.5,
    "bg_on_data": True,
    "vmax": 3,
    "vmin": -3,
    "cmap": "coolwarm",
}


def interpolate_surf_image(surf_img1, surf_img2, alpha=0.5):
    """Interpolate two surface images."""
    # Create a new surface image with the same mesh as the input
    surf_img_interpolated = deepcopy(surf_img1)
    # Interpolate the data
    for hemi in ["left", "right"]:
        surf_img_interpolated.data.parts[hemi] = (
            surf_img1.data.parts[hemi] * (1 - alpha)
            + surf_img2.data.parts[hemi] * alpha
        )
    return surf_img_interpolated


# Create figure
fig = plt.figure(figsize=(10, 8))


# Define a function to update the figure for each frame
def update(frame):
    plt.clf()

    if frame <= 10:
        # Interpolation frames (0-10)
        alpha = frame / 10
        surf_interpolated = interpolate_surf_image(
            surf_audio_source, surf_aligned, alpha=alpha
        )
        plotting.plot_surf_stat_map(
            surf_mesh=fsaverage_meshes["pial"],
            stat_map=surf_interpolated,
            figure=fig,
            **plotting_params,
        )
        plt.suptitle(
            f"Interpolated audio sentence alpha={alpha:.1f}", fontsize=16
        )
    else:
        # Target image (frame 10)
        plotting.plot_surf_stat_map(
            surf_mesh=fsaverage_meshes["pial"],
            stat_map=surf_audio_target,
            figure=fig,
            **plotting_params,
        )
        plt.suptitle("Target image", fontsize=16)

    return [fig]


# Create the animation
anim = FuncAnimation(fig, update, frames=range(12), interval=300, blit=True)

Once Loop Reflect

/home/runner/work/fmralign/fmralign/examples/plot_surf_alignment.py:182: UserWarning: Meshes are not identical but have compatible number of vertices.
  plotting.plot_surf_stat_map(
/home/runner/work/fmralign/fmralign/examples/plot_surf_alignment.py:182: DeprecationWarning: The `darkness` parameter will be deprecated in release 0.13. We recommend setting `darkness` to None
  plotting.plot_surf_stat_map(