QSIPrep

# Load QSIPrep
import module
await module.load('qsiprep/1.0.1')
await module.list()

# Import the necessary libraries
import os, base64
import json
from pathlib import Path
from IPython.display import display
from bs4 import BeautifulSoup
from IPython.display import HTML, display

# Request a freesurfer license from https://surfer.nmr.mgh.harvard.edu/registration.html
# and store it in your homedirectory 
# This is just an exampe - please replace with your license id:
!echo "Steffen.Bollmann@cai.uq.edu.au" > ~/.license
!echo "21029" >> ~/.license
!echo "*Cqyn12sqTCxo" >> ~/.license
!echo "FSxgcvGkNR59Y" >> ~/.license

%%bash
# Install the dataset from GitHub using DataLad, then fetch one subject.
if [ ! -d ds001226 ]; then
  datalad install https://github.com/OpenNeuroDatasets/ds001226.git
fi
cd ds001226 && datalad get sub-CON01

subj = Path("ds001226/sub-CON01/ses-preop")
print("=== directory layout ===")
for d in sorted(subj.iterdir()):
    print(d.name + "/")
    for f in sorted(d.iterdir()):
        print("   ", f.name)

print("\n=== DWI phase-encoding sidecars ===")
for acq in ["AP", "PA"]:
    p = subj / "dwi" / f"sub-CON01_ses-preop_acq-{acq}_dwi.json"
    j = json.loads(p.read_text())
    bval_path = subj / "dwi" / f"sub-CON01_ses-preop_acq-{acq}_dwi.bval"
    bvals = bval_path.read_text().split()
    print(f"acq-{acq}: PhaseEncodingDirection={j['PhaseEncodingDirection']!r}  "
          f"TotalReadoutTime={j['TotalReadoutTime']}  "
          f"num_volumes={len(bvals)}  bvals={bvals}")

%%bash
set -e
SUBJ_DIR="ds001226/sub-CON01/ses-preop"
mkdir -p "$SUBJ_DIR/fmap"

# Copy ONLY the .nii.gz and .json (NOT bval/bvec - not valid for _epi fieldmaps)
cp "$SUBJ_DIR/dwi/sub-CON01_ses-preop_acq-PA_dwi.nii.gz" \
   "$SUBJ_DIR/fmap/sub-CON01_ses-preop_dir-PA_epi.nii.gz"
cp "$SUBJ_DIR/dwi/sub-CON01_ses-preop_acq-PA_dwi.json" \
   "$SUBJ_DIR/fmap/sub-CON01_ses-preop_dir-PA_epi.json"

# Update the EPI sidecar: add IntendedFor, ensure PhaseEncodingDirection + TotalReadoutTime present
python3 -c "
import json
p = '$SUBJ_DIR/fmap/sub-CON01_ses-preop_dir-PA_epi.json'
j = json.load(open(p))
j['IntendedFor'] = ['ses-preop/dwi/sub-CON01_ses-preop_acq-AP_dwi.nii.gz']
assert 'PhaseEncodingDirection' in j, 'PE direction missing!'
assert 'TotalReadoutTime' in j, 'TotalReadoutTime missing!'
json.dump(j, open(p, 'w'), indent=2)
print('PE:', j['PhaseEncodingDirection'], '| TRT:', j['TotalReadoutTime'], '| IntendedFor:', j['IntendedFor'])
"

# Remove the PA from dwi/ so QSIPrep doesn't try to treat it as a DWI series to concatenate
rm "$SUBJ_DIR/dwi/sub-CON01_ses-preop_acq-PA_dwi."*

echo "Reorganisation complete. fmap/ contents:"
ls "$SUBJ_DIR/fmap/"
echo "dwi/ contents:"
ls "$SUBJ_DIR/dwi/"

%%bash
set -u -o pipefail

# FreeSurfer subjects dir must exist at startup
export SUBJECTS_DIR="$HOME/qsiprep-freesurfer-dir"
export APPTAINERENV_SUBJECTS_DIR=$SUBJECTS_DIR
export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
mkdir -p "$SUBJECTS_DIR"

# Use most (not all) of the available cores; leave 1-2 free
N_CPUS=$(nproc)
NTHREADS=$(( N_CPUS > 2 ? N_CPUS - 4 : 1 ))   # leave 1 core free
OMP_THREADS=$(( NTHREADS >= 4 ? 4 : 2 ))       # per-process cap

# Use 75% of RAM
MEM_MB=$(awk '/MemTotal/ {printf "%d", $2/1024}' /proc/meminfo)
QSIPREP_MEM=$(( MEM_MB * 75 / 100 ))

echo "Cores: ${N_CPUS} -> using ${NTHREADS} (omp ${OMP_THREADS}); RAM cap: ${QSIPREP_MEM} MB"

qsiprep ds001226 qsiprep-output participant \
    --participant-label CON01 \
    --fs-license-file ~/.license \
    --output-resolution 2 \
    -w qsiprep-work \
    --nthreads "${N_CPUS}" \
    --omp-nthreads "${OMP_THREADS}" \
    --mem-mb "${QSIPREP_MEM}" \
    2>&1 | tee qsiprep.log

%%bash
# Verify the results exists (preprocessed DWI, not just the report)
if ! ls qsiprep-output/sub-CON01/ses-preop/dwi/*desc-preproc_dwi.nii.gz >/dev/null 2>&1; then
  echo "ERROR: QSIPrep did not produce a preprocessed DWI"; exit 1
fi
echo "QSIPrep complete - preprocessed DWI found."

# Helper function to display the qsiprep report with inline figures and gifs
def display_report(report_path):

    report_dir = os.path.dirname(report_path)
    with open(report_path) as f:
        soup = BeautifulSoup(f, 'html.parser')

    # Inline <img> tags: PNG as base64, SVG, gif embedded directly
    for img in soup.find_all('img'):
        src = img.get('src', '')
        if not src:
            continue
        img_path = os.path.join(report_dir, src.lstrip('./'))
        if os.path.exists(img_path):
            if src.endswith('.svg'):                                                                                           
                with open(img_path) as f:                                                                                      
                    svg_tag = BeautifulSoup(f.read(), 'html.parser').find('svg')                                               
                    if svg_tag:                                                                                                
                        img.replace_with(svg_tag)                                                                              
            elif src.endswith('.png'):                                                                                         
                with open(img_path, 'rb') as f:                                                                                
                    img['src'] = f'data:image/png;base64,{base64.b64encode(f.read()).decode()}'                                
            elif src.endswith('.gif'):                                                                                         
                with open(img_path, 'rb') as f:                                                                                
                    img['src'] = f'data:image/gif;base64,{base64.b64encode(f.read()).decode()}'

    # Inline <object data="...svg"> tags (these were missing before!)
    for obj in soup.find_all('object'):
        data = obj.get('data', '')
        if data.endswith('.svg'):
            obj_path = os.path.join(report_dir, data.lstrip('./'))
            if os.path.exists(obj_path):
                with open(obj_path) as f:
                    svg_tag = BeautifulSoup(f.read(), 'html.parser').find('svg')
                    if svg_tag:
                        obj.replace_with(svg_tag)

    # Remove ALL tags that can execute scripts, load external resources,
    # or clutter the output. This prevents Jupyter Lab DOM interference.
    for tag in soup.find_all(['script', 'link', 'style', 'meta', 'head',
                              'nav', 'button', 'noscript', 'title']):
        tag.decompose()

    # The figures are inlined above, so drop the report's "Get figure file: <link>" lines:
    # they point into qsiprep-output/, which is not published with this page, and would 404.
    for small in soup.find_all('small'):
        if small.get_text().startswith('Get figure file'):
            small.decompose()

    # Unwrap <html> and <body> — keep their children, drop the wrapper tags
    for tag in soup.find_all(['html', 'body']):
        tag.unwrap()

    # QSIPrep's methods text links a docs page that no longer exists; point it at the
    # pipeline description in the docs for the QSIPrep version loaded above (1.0.1).
    html = str(soup).replace('https://qsiprep.readthedocs.io/en/latest/workflows.html',
                             'https://qsiprep.readthedocs.io/en/1.0.1/preprocessing.html')
    display(HTML(html))

display_report('./qsiprep-output/sub-CON01.html')

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}")