Updates a Gifti vertex' coordinates by translating them according to RAS information found in file header

translate_gifti_ras.py

# 1. Read RAS coordinates from gifti header
# and apply the underlying translation to mesh coordinates.
# 2. Save results in a new gifti files with modified header.

# %%
import nibabel as nib
import numpy as np
import os
from pathlib import Path
from tqdm import tqdm

# %%
# List gifti files to be transformed
MESH_FOLDER = Path("/home/alexis/singbrain/data/mebrains_template")
OUTPUT_FOLDER = Path(
    "/home/alexis/singbrain/data/mebrains_template_translated"
)

pial_left_path = MESH_FOLDER / "lh.MEBRAINS.pial.gii"
pial_right_path = MESH_FOLDER / "rh.MEBRAINS.pial.gii"
wm_left_path = MESH_FOLDER / "lh.MEBRAINS.smoothwm.gii"
wm_right_path = MESH_FOLDER / "rh.MEBRAINS.smoothwm.gii"

mesh_paths = [pial_left_path, pial_right_path, wm_left_path, wm_right_path]

# %%
# Check that all meshes have the same RAS transform
ras_translation = None
for mesh_path in mesh_paths:
    mesh = nib.load(mesh_path)
    translation = np.array(
        [
            float(mesh.darrays[0].meta["VolGeomC_R"]),
            float(mesh.darrays[0].meta["VolGeomC_A"]),
            float(mesh.darrays[0].meta["VolGeomC_S"]),
        ]
    )

    if ras_translation is None:
        ras_translation = translation

    assert np.all(ras_translation == translation)


# %%
def translate_mesh(mesh_path):
    mesh = nib.load(mesh_path)
    coordinates = mesh.darrays[0].data
    translation = np.array(
        [
            float(mesh.darrays[0].meta["VolGeomC_R"]),
            float(mesh.darrays[0].meta["VolGeomC_A"]),
            float(mesh.darrays[0].meta["VolGeomC_S"]),
        ]
    )

    # Compute updated coordinates
    updated_coordinates = coordinates + translation

    # Compute updated header: keep only non-RAS information
    meta = mesh.darrays[0].meta
    meta_items_to_keep = [
        "AnatomicalStructurePrimary",
        "AnatomicalStructureSecondary",
        "GeometricType",
        "Name",
    ]
    updated_meta = nib.gifti.GiftiMetaData(
        {k: meta[k] for k in meta_items_to_keep}
    )

    # Update gifti file
    mesh.darrays[0].data = updated_coordinates
    mesh.darrays[0].meta = updated_meta

    # Save updated mesh
    mesh_basename = os.path.basename(mesh_path)
    mesh.to_filename(OUTPUT_FOLDER / mesh_basename)


# %%
for mesh_path in tqdm(mesh_paths):
    translate_mesh(mesh_path)