{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Template-Space Projection\n\nWhile data is shared in a subject's own T1w space, it is\npossible to move the maps onto a shared template, making them\nready for group-level comparison.\n\nTwo routes share the same ``to_template`` entry point:\n\nThe volume route lands on MNI305 via the linear affine\nFreeSurfer fits during ``recon-all`` (stored under the historical\nname ``talairach.lta``).\n\nThe surface route walks the volume through ``vol_to_surf`` onto\nthe subject's ``fsnative`` mesh, then resamples onto\n``fsaverage``.\n\nBoth routes accept single maps, batched\n``(n_trials, n_voxels)`` arrays, and ``fsnative``-surface arrays.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>Run :doc:`plot_01 <plot_01_quickstart>` first so the shared\n   data directory has the FreeSurfer recon and anatomical mask\n   this example reads. Install the template extras with\n   ``uv sync --extra template`` or\n   ``pip install \"laion-fmri[template]\"``.</p></div>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Bind the shared data directory and load one subject\n\nPlot_01 populated the quickstart directory with the FreeSurfer\nrecon and anatomical mask that the projection chain reads.\nPick up the same directory here and load one subject.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\n\nfrom laion_fmri.config import dataset_initialize\nfrom laion_fmri.subject import load_subject\n\n# define and initialize the data directory\ndata_dir = os.environ.get(\n    \"LAION_FMRI_EXAMPLE_DATA_DIR\",\n    os.path.join(os.getcwd(), \"laion_fmri_quickstart\"),\n)\nos.makedirs(data_dir, exist_ok=True)\ndataset_initialize(data_dir)\n\n# set subject information\nsubject_id = \"sub-01\"\nsession = \"ses-01\"\nroi = \"FFA1\"\n\n# load and inspect the subject\nsub = load_subject(subject_id)\nprint(f\"Subject:        {sub.subject_id}\")\nprint(f\"FS recon ready: {sub.has_freesurfer()}\")\nprint(f\"Anatomical:     {sub.has_anatomical()}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Pick a per-voxel map to project\n\nAnything that has one value per brain-mask voxel can travel\nthrough the projection. A session-level noise-ceiling map is a\nconvenient starting point but mean betas, decoding accuracies,\nor any other per-voxel summary slot in the same way.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n\nnc = sub.get_noise_ceiling(session=session)\nprint(f\"NC shape:    {nc.shape}\")\nprint(f\"NC range:    [{np.nanmin(nc):.2f}, {np.nanmax(nc):.2f}]\")\nprint(f\"NC > 0.2:    {(nc > 0.2).sum()} voxels above threshold\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Volume route: T1w \u2192 MNI305\n\nStarting with the volume route: ``to_template`` reads the\nrecon's ``talairach.lta`` affine and resamples the\nnoise-ceiling map onto the MNI305 reference grid that\ntemplateflow ships. The return is a 3-D NIfTI, ready for\nfurther handling. Here, the focus is on visualization.\nVolumetric alignment from a single linear affine is good\nenough for whole-brain visualizations, but cortical analyses\nthat need sub-millimetre accuracy will do better with the\nsurface route further down.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mni305_img = sub.to_template(nc, \"MNI305\")\nprint(f\"MNI305 image: shape={mni305_img.shape}, \"\n      f\"dtype={mni305_img.get_data_dtype()}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualize the MNI305 output\n\nA good way to convince yourself the projection worked is to\noverlay the resampled map on the matching MNI305 anatomy and\nsee whether high noise-ceiling voxels land where they should.\nBuilding that figure also makes a nice excuse to walk through\nhow ``laion_fmri`` plays with the wider neuroimaging\necosystem, so the code below is intentionally a little more\nverbose than strictly necessary. The plan is to (i) pull the\nMNI305 T1w and brain mask from ``templateflow``, (ii) use\n``nilearn`` to multiply the T1w by the brain mask and crop\nto the brain's bounding box so the backdrop shows cortex\nrather than skull and air, (iii) threshold the noise-ceiling\nmap at 10% variance explained, and (iv) overlay the result on\nthe prepared anatomy. This is the typical pattern for any\ngroup-level figure on this dataset.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport seaborn as sns\nfrom matplotlib.colors import Normalize\nfrom nilearn import plotting\nfrom nilearn.image import crop_img, math_img\nfrom templateflow.api import get as tflow_get\n\n# use templateflow to get the mni305 template\nmni305_t1w = tflow_get(\n    \"MNI305\", suffix=\"T1w\", extension=\".nii.gz\",\n)\nmni305_brain_mask = tflow_get(\n    \"MNI305\", suffix=\"mask\", desc=\"brain\", extension=\".nii.gz\",\n)\n\n# Templateflow ships a head T1w and a brain mask for MNI305 but\n# no brain-extracted T1w directly. Thus, the T1w is multiplied\n# by the brain mask to drop the skull / dura / air and then\n# cropped to the brain's bounding box so plot_stat_map renders\n# only the cortex.\nmni305_bg = crop_img(\n    math_img(\n        \"img * mask\",\n        img=str(mni305_t1w), mask=str(mni305_brain_mask),\n    )\n)\n\n# define the colormap and data range\nmako_cmap = sns.color_palette(\"mako\", as_cmap=True)\nnc_vmax = float(nc.max())\n\n# define the NC threshold and cuts to plot. The threshold is\n# shared with the surface figure further down so colors read\n# the same way on both.\nthreshold = 10.0\ncuts = [-17, -5, 8]\n\n# set up the figure\nfig = plt.figure(figsize=(10, 4), facecolor=\"white\")\ngs = fig.add_gridspec(2, 1, height_ratios=[1, 0.05], hspace=0.1)\nfig.subplots_adjust(top=0.98, bottom=0.15)\nstrip_gs = gs[0].subgridspec(1, 3, wspace=0.05)\naxes = [fig.add_subplot(strip_gs[0, i]) for i in range(3)]\n\n# set up the colorbar strip, narrowed to match the surface\n# figure's colorbar width-vs-figure-width proportion\ncbar_strip_gs = gs[1].subgridspec(\n    1, 3, width_ratios=[0.15, 0.7, 0.15],\n)\ncbar_ax = fig.add_subplot(cbar_strip_gs[0, 1])\n\n# plot the different cuts\nfor ax, z in zip(axes, cuts):\n    ax.set_facecolor(\"white\")\n    plotting.plot_stat_map(\n        mni305_img, bg_img=mni305_bg, axes=ax,\n        display_mode=\"z\", cut_coords=[z],\n        cmap=mako_cmap, vmax=nc_vmax, colorbar=False,\n        black_bg=False, threshold=threshold,\n    )\n\n# define and render the colorbar\nsm = plt.cm.ScalarMappable(\n    cmap=mako_cmap, norm=Normalize(vmin=0, vmax=nc_vmax),\n)\nfig.colorbar(\n    sm, cax=cbar_ax, orientation=\"horizontal\",\n    label=f\"{session} noise ceiling on MNI305 (% var. expl.)\",\n)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Surface route: T1w \u2192 fsaverage (single hemisphere)\n\nSurface projections take a two-step trip: ``vol_to_surf`` lifts\nthe volume onto the subject's ``fsnative`` mesh, then a surface\nresampler carries it across to ``fsaverage``. The\n``hemisphere`` can be picked with ``hemi=\"L\"`` or ``\"R\"`` and\nthe ``mesh density`` with ``fsaverage_density``. The default\n``fsaverage5`` (10k vertices per hemi) is fast; ``fsaverage6``\n(41k) and ``fsaverage`` (164k) trade compute time for finer\ndetail.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nc_fsavg_lh = sub.to_template(nc, \"fsaverage\", hemi=\"L\")\nnc_fsavg_rh = sub.to_template(nc, \"fsaverage\", hemi=\"R\")\nprint(f\"fsaverage5 L: {nc_fsavg_lh.shape}\")\nprint(f\"fsaverage5 R: {nc_fsavg_rh.shape}\")\n\n# Higher density (fsaverage6 / 41k vertices per hemi)\nnc_fsavg6_lh = sub.to_template(\n    nc, \"fsaverage\", hemi=\"L\", fsaverage_density=\"fsaverage6\",\n)\nprint(f\"fsaverage6 L: {nc_fsavg6_lh.shape}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Re-project the noise ceiling at the full fsaverage density\n\nTo showcase the high resolution in the visualization, the\nprojection is re-run to ``fsaverage`` (164k vertices per hemi).\nHowever, the choice of the surface depends on the intended use\ncase and analyses.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nc_fsavg_lh_hi = sub.to_template(\n    nc, \"fsaverage\", hemi=\"L\", fsaverage_density=\"fsaverage\",\n)\nnc_fsavg_rh_hi = sub.to_template(\n    nc, \"fsaverage\", hemi=\"R\", fsaverage_density=\"fsaverage\",\n)\nprint(f\"fsaverage L (164k): {nc_fsavg_lh_hi.shape}\")\nprint(f\"fsaverage R (164k): {nc_fsavg_rh_hi.shape}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Configure the gyri / sulci backdrop and the view layout\n\nTwo things make a surface figure read well: a clear anatomical\nbackdrop (so the reader can tell where on the cortex they are\nlooking) and a layout that shows enough viewpoints to cover\nthe regions of interest. The cell below sets up both at once.\nThe backdrop encodes gyri / sulci in two shades of gray so\nthe noise-ceiling overlay stays the most colorful element on\nthe plot, and the view layout is a per-row tuple that pairs\neach anatomical view (``lateral``, ``medial``, ``posterior``,\n``flat``) with the matching ``fsaverage`` mesh. This is the\nsame recipe to reach for whenever a surface figure has to be\ncompared cleanly against a volume one further down.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from matplotlib.colors import ListedColormap, to_rgba\nfrom nilearn import datasets, surface\n\n# define colors for the surface plots\nGYRI_HEX = \"#9B978D\"\nSULCI_HEX = \"#595959\"\nN_NC_BINS = 64\nfsavg = datasets.fetch_surf_fsaverage(\"fsaverage\")  # 164k\n\n# stack the gyri / sulci greys ahead of 64 mako bins so a single\n# integer label per vertex covers both anatomy and data\nmako_bins = [\n    mako_cmap(i / (N_NC_BINS - 1)) for i in range(N_NC_BINS)\n]\ncomposite_cmap = ListedColormap(\n    [(1, 1, 1, 1), to_rgba(GYRI_HEX), to_rgba(SULCI_HEX)]\n    + mako_bins\n)\n\n# define a per-row layout with multiple views:\n# (label, (lh_view, rh_view), lh_mesh, rh_mesh).\nview_rows = [\n    (\"lateral\",\n     (\"lateral\", \"lateral\"),\n     fsavg.infl_left, fsavg.infl_right),\n    (\"medial\",\n     (\"medial\", \"medial\"),\n     fsavg.infl_left, fsavg.infl_right),\n    (\"posterior\",\n     (\"posterior\", \"posterior\"),\n     fsavg.infl_left, fsavg.infl_right),\n    (\"flat\",\n     ((90, -90), (90, -90)),\n     fsavg.flat_left, fsavg.flat_right),\n]\n\n# define a shared spatial extent for both flat panels.\n# Without this, matplotlib 3D auto-fits each panel to its own\n# mesh and the two end up at slightly different vertical\n# positions\nlh_flat_verts = surface.load_surf_mesh(fsavg.flat_left)[0]\nrh_flat_verts = surface.load_surf_mesh(fsavg.flat_right)[0]\nflat_all = np.concatenate([lh_flat_verts, rh_flat_verts])\nflat_xlim = (flat_all[:, 0].min(), flat_all[:, 0].max())\nflat_ylim = (flat_all[:, 1].min(), flat_all[:, 1].max())\nflat_zlim = (flat_all[:, 2].min(), flat_all[:, 2].max())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the multi-view surface grid\n\nWith the layout and the composite colormap in place, the goal\nof this cell is to render four anatomical views per hemisphere\n(lateral, medial, posterior, flat) using the *same* noise-\nceiling threshold and color range as the volume figure above.\nKeeping both figures on the same scale matters because it lets\nthe reader make a direct visual comparison: a noise-ceiling\nvalue of, say, 0.4 produces the exact same color whether it\nis shown on the MNI305 slices or on the inflated fsaverage\nsurface, and the thresholded boundary marks the same cut on\nboth sides. Without that alignment, two figures of the same\nunderlying data can look misleadingly different.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# setup the figure\nfig = plt.figure(figsize=(9.0, 11.5), facecolor=\"white\")\ngs = fig.add_gridspec(\n    5, 1, height_ratios=[1, 1, 1, 1, 0.06], hspace=0.02,\n)\nfig.subplots_adjust(top=0.99, bottom=0.07, left=0.04, right=0.99)\n\n# loop over views\nfor row_idx, (label, views, mesh_l, mesh_r) in enumerate(view_rows):\n    row_gs = gs[row_idx].subgridspec(1, 2, wspace=0.02)\n    for col_idx, (mesh, arr, curv_path, hemi_full, view) in (\n        enumerate(zip(\n            (mesh_l, mesh_r),\n            (nc_fsavg_lh_hi, nc_fsavg_rh_hi),\n            (fsavg.curv_left, fsavg.curv_right),\n            (\"left\", \"right\"),\n            views,\n        ))\n    ):\n        ax = fig.add_subplot(row_gs[0, col_idx], projection=\"3d\")\n        curv = surface.load_surf_data(curv_path)\n        roi_map = np.where(curv < 0, 1.0, 2.0)\n        above = arr >= threshold\n        if above.any():\n            norm = np.clip(arr / nc_vmax, 0.0, 1.0)\n            bin_idx = (norm * (N_NC_BINS - 1)).astype(np.int32)\n            roi_map[above] = 3.0 + bin_idx[above]\n        plotting.plot_surf_roi(\n            mesh, roi_map=roi_map,\n            hemi=hemi_full, view=view,\n            cmap=composite_cmap, vmin=0, vmax=2 + N_NC_BINS,\n            colorbar=False, figure=fig, axes=ax,\n        )\n        ax.set_anchor(\"C\")\n        if label == \"flat\":\n            ax.set_xlim(flat_xlim)\n            ax.set_ylim(flat_ylim)\n            ax.set_zlim(flat_zlim)\n        if col_idx == 0:\n            ax.text2D(\n                -0.05, 0.5, label,\n                transform=ax.transAxes,\n                rotation=90, ha=\"center\", va=\"center\",\n                fontsize=11,\n            )\n\n# define and render the colorbar\ncbar_strip_gs = gs[4].subgridspec(\n    1, 3, width_ratios=[0.15, 0.7, 0.15],\n)\ncbar_ax = fig.add_subplot(cbar_strip_gs[0, 1])\nsm_surf = plt.cm.ScalarMappable(\n    cmap=mako_cmap,\n    norm=Normalize(vmin=0, vmax=nc_vmax),\n)\nfig.colorbar(\n    sm_surf, cax=cbar_ax, orientation=\"horizontal\",\n    label=f\"{session} noise ceiling on fsaverage (% var. expl.)\",\n)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Both hemispheres at once (``hemi=None``)\n\nIf both hemispheres should be returned in a single call,\n``hemi`` can be left at its default of ``None``. The return\nvalue becomes a dict keyed by ``\"L\"`` and ``\"R\"``, with each\nvalue the same 1-D array shape as the single-hemi call above.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# project both hemispheres at once\nboth_hemis = sub.to_template(nc, \"fsaverage\")\nprint(f\"Returned: {type(both_hemis).__name__} with keys \"\n      f\"{sorted(both_hemis)}\")\nprint(f\"L shape: {both_hemis['L'].shape}\")\nprint(f\"R shape: {both_hemis['R'].shape}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Batched input: (n_trials, n_voxels) \u2192 (n_trials, n_vertices)\n\nSingle-trial analyses usually start from a ``(n_trials,\nn_voxels)`` array rather than a single summary map. The same\n``to_template`` call accepts the 2-D input and keeps the trial\naxis on the output: a 4-D NIfTI for volume targets, a 2-D\n``(n_trials, n_vertices)`` array for surface targets. Inputs\nare sized to the brain mask. For ROI-masked betas, scatter\nthem back into the brain-mask shape first with\n:meth:`Subject.to_nifti`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# load the first eight single-trial betas at the brain-mask\n# resolution. ``streaming=True`` keeps peak memory low while\n# the session NIfTI is decompressed and read volume by volume.\nbetas = sub.get_betas(session=session, streaming=True)\nbatch_full = betas[:2]\nprint(f\"Input batch:  {batch_full.shape}\")\n\n# project the batched input onto MNI305 (4-D NIfTI, trial axis\n# trailing)\nmni305_batch = sub.to_template(batch_full, \"MNI305\")\nprint(f\"MNI305 4-D:   shape={mni305_batch.shape}\")\n\n# project the batched input onto fsaverage5 (2-D, trial axis\n# leading)\nfsavg_batch = sub.to_template(batch_full, \"fsaverage\", hemi=\"L\")\nprint(f\"fsavg5 batch: {fsavg_batch.shape}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Surface \u2192 surface (fsnative \u2192 fsaverage)\n\nIf the data already lives on the subject's ``fsnative`` mesh\n(e.g., a surface-based ROI mask, surface-fit betas, or\nanything else in surface space), ``surface_to_template``\nskips the volume step and resamples straight onto\n``fsaverage``. Pass a per-hemi array with ``hemi=\"L\"`` or\n``\"R\"``, or a ``{\"L\": ..., \"R\": ...}`` dict for both\nhemispheres in one call. The input vertex count has to match\nthe subject's ``fsnative`` mesh.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# pull an fsnative-space ROI mask from disk\nroi_data = sub.get_roi_data(roi, format=\"func.gii\", hemi=\"L\")\nfsnative_lh = roi_data[roi][\"gii\"][\"hemi-L\"][\"func.gii\"]\nprint(f\"fsnative L mesh size: {fsnative_lh.shape}\")\n\n# resample the fsnative ROI mask onto fsaverage5\nffa1_on_fsavg = sub.surface_to_template(\n    fsnative_lh.astype(np.float32), hemi=\"L\",\n)\nprint(f\"FFA1 on fsavg5 L:     {ffa1_on_fsavg.shape}\")\nprint(f\"Vertices marked ROI:  {int((ffa1_on_fsavg > 0.5).sum())}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Direction-specific entry points\n\n``to_template`` infers the route from the target name. If the\ninput / output direction should be visible at the call site,\nthree direction-specific methods do the same work under more\ndescriptive names:\n\n* ``volume_to_template(values, target)``: volume in, volume\n  out (e.g. T1w \u2192 MNI305).\n* ``volume_to_surface(values, target=\"fsaverage\")``: volume\n  in, surface out (e.g. T1w \u2192 fsaverage).\n* ``surface_to_template(values, target=\"fsaverage\", hemi=...)``:\n  surface in, surface out (e.g. fsnative \u2192 fsaverage).\n\nA target that doesn't match the chosen direction raises\n``ValueError`` at the call site, so a wrong combination is\ncaught before any data moves.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# direct volume-to-volume projection\nvol_mni = sub.volume_to_template(nc, \"MNI305\")\n# direct volume-to-surface projection\nvol_surf_l = sub.volume_to_surface(nc, hemi=\"L\")\nprint(f\"volume_to_template(MNI305):   {vol_mni.shape}\")\nprint(f\"volume_to_surface(L):         {vol_surf_l.shape}\")\n\n# For example, when the target doesn't match the chosen\n# direction (uncomment locally to see the raised ValueError;\n# kept commented so the gallery build does not crash):\n# sub.volume_to_template(nc, \"fsaverage\")  # ValueError"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}