[WIP] add plot_anat_landmarks function

doc/Makefile

@@ -9,6 +9,11 @@ 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)
@@ -23,3 +28,11 @@ help:
 # View the built documentation
 view:
 	@python -c "import webbrowser; webbrowser.open_new_tab('file://$(PWD)/_build/html/index.html')"
+
+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

@@ -53,7 +53,6 @@ mne_bids.stats
 
    count_events
 
-
 mne_bids.copyfiles
 ------------------
 
@@ -74,3 +73,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_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
@@ -134,16 +135,20 @@
 anat_dir = t1w_bids_path.directory
 
 ###############################################################################
-# Let's have another look at our BIDS directory
+# 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)
+plt.suptitle('T1 MRI')
+
 ###############################################################################
 # Our BIDS dataset is now ready to be shared. We can easily estimate the
-# transformation matrix using ``MNE-BIDS`` and the BIDS dataset.
+# 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)
 
 ###############################################################################
-# Finally, let's use the T1 weighted MRI image and plot the anatomical
+# One can also directly use the T1 weighted MRI image and plot the anatomical
 # landmarks Nasion, LPA, and RPA onto the brain image. For that, we can
 # extract the location of Nasion, LPA, and RPA from the MEG file, apply our
 # transformation matrix :code:`trans`, and plot the results.
@@ -169,7 +174,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, :],
@@ -214,9 +219,8 @@
 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)
+plt.suptitle('T1 Defaced')
 
 ###############################################################################
 # Writing defaced and anonymized FLASH MRI image
@@ -250,9 +254,8 @@
 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)
+plt.suptitle('Flash Defaced')
 
 ###############################################################################
 # Option 2 : Use manual landmarks coordinates in scanner RAS for FLASH image
@@ -282,9 +285,8 @@
 )
 
 # 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)
+plt.suptitle('Flash Defaced')
 
 ###############################################################################
 # Option 3 : Compute the landmarks in scanner RAS or mri voxel space from trans
@@ -324,9 +326,8 @@
 )
 
 # 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)
+plt.suptitle('Flash Defaced')
 
 ##############################################################################
 # Let's now pass it in voxel coordinates
@@ -352,9 +353,8 @@
 )
 
 # 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)
+plt.suptitle('Flash Defaced')
 
 ###############################################################################
 # .. LINKS

mne_bids/__init__.py

@@ -12,3 +12,4 @@
                             write_meg_calibration, write_meg_crosstalk)
 from mne_bids.sidecar_updates import update_sidecar_json
 from mne_bids.inspect import inspect_dataset
+from mne_bids import viz

mne_bids/tests/test_viz.py

@@ -0,0 +1,36 @@
+from pathlib import Path
+
+import pytest
+
+import mne
+from mne.datasets import testing
+from mne_bids.viz import plot_anat_landmarks
+from mne_bids import BIDSPath, write_anat
+
+from mne.utils import requires_nibabel, requires_version
+
+
+@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)
+
+    bids_root = Path(tmpdir)
+    t1w_mgh_fname = Path(data_path) / 'subjects' / '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.warns(RuntimeWarning, match='No landmarks available'):
+        fig = plot_anat_landmarks(bids_path, show=False)
+        assert fig is None
+
+    bids_path = write_anat(t1w_mgh_fname, bids_path=bids_path, raw=raw,
+                           trans=trans, deface=True, verbose=True,
+                           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,70 @@
+import json
+
+import numpy as np
+from scipy import linalg
+
+import nibabel as nib
+from mne.utils import warn
+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
+        Call pyplot.show() at the end. Defaults to True.
+
+    Returns
+    -------
+    fig : matplotlib.figure.Figure | None
+        The figure object containing the plot. None if no landmarks
+        are avaiable.
+    """
+    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:
+        warn("No landmarks available with the image")
+        return
+
+    # Move the coords_dict from MRI Voxel to RAS
+    mgh = nib.MGHImage(nii.dataobj, nii.affine)
+    ras2vox = mgh.header.get_ras2vox()
+    vox2ras = linalg.inv(ras2vox)
+    for label in coords_dict:
+        vox_pos = np.array(coords_dict[label])
+        ras_pos = mne.transforms.apply_trans(vox2ras, 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,
+        )
+
+    if show:
+        plt.show()
+
+    return fig