[WIP] add plot_anat_landmarks function

doc/Makefile

@@ -9,11 +9,16 @@ PAPER         =
 SOURCEDIR     = .
 BUILDDIR      = _build
 
+# Internal variables.
+PAPEROPT_a4     = -D latex_paper_size=a4
+PAPEROPT_letter = -D latex_paper_size=letter
+ALLSPHINXOPTS   = -d _build/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) .
+
 # Put it first so that "make" without argument is like "make help".
 help:
 	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
 
-.PHONY: help Makefile view
+.PHONY: help Makefile view clean html-noplot
 
 # Catch-all target: route all unknown targets to Sphinx using the new
 # "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
@@ -26,3 +31,8 @@ view:
 
 clean:
 	-rm -rf _build auto_examples generated
+
+html-noplot:
+	@$(SPHINXBUILD) -D plot_gallery=0 -b html $(ALLSPHINXOPTS) _build/html
+	@echo
+	@echo "Build finished. The HTML pages are in _build/html."

doc/api.rst

@@ -61,7 +61,6 @@ mne_bids.stats
 
    count_events
 
-
 mne_bids.copyfiles
 ------------------
 
@@ -82,3 +81,19 @@ mne_bids.copyfiles
    copyfile_ctf
    copyfile_bti
    copyfile_kit
+
+mne_bids.viz
+--------------
+
+:py:mod:`mne_bids.viz`:
+
+.. automodule:: mne_bids.viz
+   :no-members:
+   :no-inherited-members:
+
+.. currentmodule:: mne_bids.viz
+
+.. autosummary::
+   :toctree: generated/
+
+   plot_anat_landmarks

examples/convert_mri_and_trans.py

@@ -51,6 +51,7 @@
 
 from mne_bids import (write_raw_bids, BIDSPath, write_anat, get_anat_landmarks,
                       get_head_mri_trans, print_dir_tree)
+from mne_bids.viz import plot_anat_landmarks
 
 # %%
 # We will be using the `MNE sample data <mne_sample_data_>`_ and write a basic
@@ -144,9 +145,18 @@
 anat_dir = t1w_bids_path.directory
 
 # %%
-# Let's have another look at our BIDS directory
+# Our BIDS dataset is now ready to be shared. We can easily estimate the
+# transformation matrix using ``MNE-BIDS`` and the BIDS dataset since we
+# have now anatomical landmarks.
+estim_trans = get_head_mri_trans(
+    bids_path=bids_path, fs_subject='sample', fs_subjects_dir=fs_subjects_dir)
+
+# %%
+# Let's have another look at our BIDS directory and plot the written landmarks
 print_dir_tree(output_path)
 
+plot_anat_landmarks(t1w_bids_path, vmax=160)
+
 # %%
 # Our BIDS dataset is now ready to be shared. We can easily estimate the
 # transformation matrix using ``MNE-BIDS`` and the BIDS dataset.
@@ -189,7 +199,7 @@
 t1_nii_fname = op.join(anat_dir, 'sub-01_ses-01_T1w.nii.gz')
 
 # Plot it
-fig, axs = plt.subplots(3, 1, figsize=(7, 7), facecolor='k')
+fig, axs = plt.subplots(3, 1, figsize=(7, 7))
 for point_idx, label in enumerate(('LPA', 'NAS', 'RPA')):
     plot_anat(t1_nii_fname, axes=axs[point_idx],
               cut_coords=mri_pos[point_idx, :],
@@ -232,9 +242,7 @@
 t1_nii_fname = op.join(anat_dir, 'sub-01_ses-01_T1w.nii.gz')
 
 # Plot it
-fig, ax = plt.subplots()
-plot_anat(t1_nii_fname, axes=ax, title='Defaced', vmax=160)
-plt.show()
+plot_anat_landmarks(t1w_bids_path, vmax=160)
 
 # %%
 # Writing defaced and anonymized FLASH MRI image
@@ -265,9 +273,7 @@
 flash_nii_fname = op.join(anat_dir, 'sub-01_ses-01_FLASH.nii.gz')
 
 # Plot it
-fig, ax = plt.subplots()
-plot_anat(flash_nii_fname, axes=ax, title='Defaced', vmax=700)
-plt.show()
+plot_anat_landmarks(flash_bids_path, vmax=700)
 
 # %%
 # Using manual landmark coordinates in scanner RAS
@@ -302,9 +308,7 @@
 )
 
 # Plot it
-fig, ax = plt.subplots()
-plot_anat(flash_nii_fname, axes=ax, title='Defaced', vmax=700)
-plt.show()
+plot_anat_landmarks(flash_bids_path, vmax=700)
 
 # %%
 # .. LINKS

mne_bids/__init__.py

@@ -14,5 +14,6 @@
                             get_anat_landmarks, anonymize_dataset)
 from mne_bids.sidecar_updates import update_sidecar_json, update_anat_landmarks
 from mne_bids.inspect import inspect_dataset
+from mne_bids import viz
 from mne_bids.dig import (template_to_head, convert_montage_to_ras,
                           convert_montage_to_mri)

mne_bids/tests/test_viz.py

@@ -0,0 +1,38 @@
+from pathlib import Path
+
+import pytest
+
+import mne
+from mne.datasets import testing
+from mne.utils import requires_nibabel, requires_version
+from mne_bids.viz import plot_anat_landmarks
+from mne_bids import BIDSPath, write_anat, get_anat_landmarks
+
+
+@requires_nibabel()
+@requires_version('nilearn', '0.6')
+def test_plot_anat_landmarks(tmpdir):
+    """Test writing anatomical data with pathlib.Paths."""
+    data_path = Path(testing.data_path())
+    raw_fname = data_path / 'MEG' / 'sample' / 'sample_audvis_trunc_raw.fif'
+    trans_fname = str(raw_fname).replace('_raw.fif', '-trans.fif')
+    raw = mne.io.read_raw_fif(raw_fname)
+    trans = mne.read_trans(trans_fname)
+    fs_subjects_dir = Path(data_path) / 'subjects'
+
+    bids_root = Path(tmpdir)
+    t1w_mgh_fname = fs_subjects_dir / 'sample' / 'mri' / 'T1.mgz'
+    bids_path = BIDSPath(subject='01', session='mri', root=bids_root)
+    bids_path = write_anat(t1w_mgh_fname, bids_path=bids_path, overwrite=True)
+
+    with pytest.raises(ValueError, match='No landmarks available'):
+        plot_anat_landmarks(bids_path, show=False)
+
+    landmarks = get_anat_landmarks(
+        t1w_mgh_fname, raw.info, trans, fs_subject='sample',
+        fs_subjects_dir=fs_subjects_dir)
+
+    bids_path = write_anat(t1w_mgh_fname, bids_path=bids_path,
+                           landmarks=landmarks, overwrite=True)
+    fig = plot_anat_landmarks(bids_path, show=False)
+    assert len(fig.axes) == 12  # 3 subplots + 3 x 3 MRI slices

mne_bids/viz.py

@@ -0,0 +1,64 @@
+import json
+
+import numpy as np
+import nibabel as nib
+
+import mne
+
+
+def plot_anat_landmarks(bids_path, vmax=None, show=True):
+    """Plot anatomical landmarks attached to an MRI image.
+
+    Parameters
+    ----------
+    bids_path : mne_bids.BIDSPath
+        Path of the MRI image.
+    vmax : float
+        Maximum colormap value.
+    show : bool
+        Whether to show the figure after plotting. Defaults to ``True``.
+
+    Returns
+    -------
+    fig : matplotlib.figure.Figure
+        The figure object containing the plot.
+    """
+    import matplotlib.pyplot as plt
+    from nilearn.plotting import plot_anat
+
+    nii = nib.load(str(bids_path))
+
+    json_path = bids_path.copy().update(extension=".json")
+
+    n_landmarks = 0
+    if json_path.fpath.exists():
+        json_content = json.load(open(json_path))
+        coords_dict = json_content.get("AnatomicalLandmarkCoordinates", dict())
+        n_landmarks = len(coords_dict)
+
+    if not n_landmarks:
+        raise ValueError("No landmarks available with the image")
+
+    for label in coords_dict:
+        vox_pos = np.array(coords_dict[label])
+        ras_pos = mne.transforms.apply_trans(nii.affine, vox_pos)
+        coords_dict[label] = ras_pos
+
+    ########################################################################
+    # Plot it with nilearn
+    fig, axs = plt.subplots(
+        n_landmarks, 1, figsize=(6, 2.3 * n_landmarks),
+        facecolor="w")
+
+    for point_idx, (label, ras_pos) in enumerate(coords_dict.items()):
+        plot_anat(
+            str(bids_path), axes=axs[point_idx], cut_coords=ras_pos,
+            title=label, vmax=vmax,
+        )
+
+    plt.suptitle(bids_path.fpath.name)
+
+    if show:
+        fig.show()
+
+    return fig