FreeSurfer: recon-all-clinical

# we can use module to load freesurfer in a specific version
import module
await module.load('freesurfer/8.0.0')
await module.list()

%%capture
!pip install nibabel numpy

!recon-all-clinical.sh --help

! echo "Downloading File"; mkdir -p ./data/structural && osf -p bt4ez fetch osfstorage/TOMCAT_DIB/sub-01/ses-01_7T/anat/sub-01_ses-01_7T_T1w_defaced.nii.gz ./data/structural/mp2rage.nii.gz

!ls ./data/structural/mp*

# Display original file in grayscale
from ipyniivue import NiiVue
from IPython.display import display, Image
import nibabel as nib, numpy as np

T1_PATH = "./data/structural/mp2rage.nii.gz" 

# Robust intensity window (ignore zeros and outliers)
t1 = nib.load(T1_PATH).get_fdata()
vals = t1[t1 > 0]
vmin, vmax = np.percentile(vals, [0.5, 99.5]) if vals.size else (None, None)

nv = NiiVue(
    height=600,
    multiplanar_layout="GRID",
    multiplanar_show_render=True,
    is_ruler=False,
    is_colorbar=True,
    is_orient_cube=True,
    is_radiological_convention=True,
    back_color=(0.1,0.1,0.1,1.0),
)

nv.add_volume({
    "path": T1_PATH,
    "name": "native",
    "opacity": 1.0,
    "colormap": "gray",
    "cal_min": float(vmin) if vmin is not None else None,
    "cal_max": float(vmax) if vmax is not None else None,
})

nv.set_crosshair_color([0,0,0,1])
display(nv)

%%bash
#--------------------------------------------
# FreeSurfer recon-all-clinical pipeline
#--------------------------------------------

set -euo pipefail  # Stop on unset variables or errors

#----------------------------
# Config
#----------------------------
INPUT_SCAN=./data/structural/mp2rage.nii.gz
SUBJECT_ID=TestSubject
THREADS=2 #adapt if more resources are available
SUBJECT_DIR=./data/structural/recon-all-clinical

# Ensure output folder exists
mkdir -p "$SUBJECT_DIR"
echo "SUBJECT_DIR=$SUBJECT_DIR"
echo "SUBJECT_ID=$SUBJECT_ID"
echo "INPUT_SCAN=$INPUT_SCAN"
echo

#----------------------------
# Allow deeper folder structures in some FreeSurfer builds
#----------------------------
export FS_ALLOW_DEEP=1
export SINGULARITYENV_FS_ALLOW_DEEP=1
export APPTAINERENV_FS_ALLOW_DEEP=1

#----------------------------
# Start pipeline
#----------------------------
# Remove previous "is-running" file if it exists
rm -f "$SUBJECT_DIR/$SUBJECT_ID/scripts/IsRunning.lh+rh"

# Run the actual recon-all-clinical pipeline
recon-all-clinical.sh "$INPUT_SCAN" "$SUBJECT_ID" "$THREADS" "$SUBJECT_DIR"

#----------------------------
# End message
#----------------------------
echo "----- recon-all-clinical.sh completed -----"
echo "Results saved to: $SUBJECT_DIR/$SUBJECT_ID"

# Display aseg.mgz in FreeSurfer colors, labels 0–100
import tempfile, os
from pathlib import Path

ASEG_PATH = "./data/structural/recon-all-clinical/TestSubject/mri/aseg.mgz"
assert Path(ASEG_PATH).exists(), f"Missing aseg: {ASEG_PATH}"

aseg_img = nib.load(ASEG_PATH)
aseg = np.rint(aseg_img.get_fdata()).astype(np.int32)

# Keep only labels 0–100, everything else to 0
aseg = np.where((aseg >= 0) & (aseg <= 100), aseg, 0)

# Save temporary integer NIfTI
tmp_path = os.path.join(tempfile.gettempdir(), "aseg_fs_0to100.nii.gz")
nib.save(nib.Nifti1Image(aseg, aseg_img.affine, aseg_img.header), tmp_path)

# Viewer
nv = NiiVue(
    height=600,
    multiplanar_layout="GRID",
    multiplanar_show_render=True,
    is_colorbar=True,
    is_orient_cube=True,
    is_radiological_convention=True,
    back_color=(0.1, 0.1, 0.1, 1.0),
)

# Use FreeSurfer colormap
nv.add_volume({
    "path": tmp_path,
    "name": "aseg",
    "opacity": 1.0,
    "colormap": "freesurfer",
})

nv.set_crosshair_color([0,0,0,1])
display(nv)

# Pial surfaces with explicit per-mesh color + shader (ipyniivue)
from ipyniivue import Mesh

SUBJECTS_DIR = "./data/structural/recon-all-clinical"
SUBJECT = "TestSubject"
surf_dir = Path(SUBJECTS_DIR) / SUBJECT / "surf"

lh = surf_dir / "lh.pial"
rh = surf_dir / "rh.pial"
assert lh.exists() and rh.exists(), "Missing pial surfaces"

nv = NiiVue(
    height=600,
    is_colorbar=False,
    is_orient_cube=True,
    is_radiological_convention=True,
    back_color=(0.05, 0.05, 0.05, 1.0),
)

# Use rgba255 (0..255 ints) — more reliable across builds than 0..1 floats
nv.add_mesh(Mesh(path=str(lh), name="lh.pial", rgba255=[255, 90, 90, 255], opacity=1.0))
nv.add_mesh(Mesh(path=str(rh), name="rh.pial", rgba255=[90, 90, 255, 255], opacity=1.0))


display(nv)

# White surfaces (lh/rh) with solid colors

lh = surf_dir/"lh.white"; rh = surf_dir/"rh.white"
assert lh.exists() and rh.exists(), "Missing white surfaces"

nv = NiiVue(height=600,is_colorbar=False, is_orient_cube=True, 
            is_radiological_convention=True, back_color=(0.05,0.05,0.05,1.0))

nv.add_mesh(Mesh(path=str(lh), name="lh.white", rgba255=[255,200,90,255], opacity=1.0))
nv.add_mesh(Mesh(path=str(rh), name="rh.white", rgba255=[90,200,255,255],  opacity=1.0))

display(nv)

# Inflated surfaces (lh/rh) with solid colors

lh = surf_dir/"lh.inflated"; rh = surf_dir/"rh.inflated"
assert lh.exists() and rh.exists(), "Missing inflated surfaces"

nv = NiiVue(height=600, is_colorbar=False, is_orient_cube=True, 
            is_radiological_convention=True, back_color=(0.05,0.05,0.05,1.0))

nv.add_mesh(Mesh(path=str(lh), name="lh.inflated", rgba255=[255,140,140,255], opacity=1.0))
#nv.add_mesh(Mesh(path=str(rh), name="rh.inflated", rgba255=[140,140,255,255], opacity=1.0))

display(nv)

# Show directory layout and recon logs
import os, subprocess, textwrap

SUBJECTS_DIR = "./data/structural/recon-all-clinical"
SUBJECT = "TestSubject"
base = f"{SUBJECTS_DIR}/{SUBJECT}"

print("=== Subject directory tree (top) ===")
subprocess.run(["bash","-lc", f'ls -lh {base}'])
print("\n=== Logs/scripts ===")
scripts = f"{base}/scripts"
subprocess.run(["bash","-lc", f'ls -lh {scripts} 2>/dev/null || true'])

notes = [
    ("mri/", "All volumetric outputs (.mgz)."),
    ("surf/", "Cortical meshes and annotation files (.white, .pial, .inflated, .annot)."),
    ("label/", "Labels, annotation tables (.ctab), region lists."),
    ("stats/", "Cortex thickness/surface area stats per region (for group analysis)."),
    ("scripts/", "Run logs, IsRunning flags, and command histories."),
]
for d, t in notes:
    print(f"{d:12s} {t}")

import os

%load_ext watermark

%watermark
%watermark --iversions

neurodesktop_version = (
    os.environ.get('JUPYTER_IMAGE', '').split(':')[-1] or
    os.environ.get('NEURODESKTOP_VERSION', 'unknown')
)

print(f"Neurodesktop version: {neurodesktop_version}")