[ENH, MRG] Use the new Dipy transform_points function for warping intracranial contacts (mne-tools#11572)

Author: alexrockhill

doc/changes/0.24.inc

@@ -137,7 +137,7 @@ Enhancements
 
 - The :meth:`mne.Epochs.get_data` method now has a ``units`` parameter (:gh:`9553` by `Stefan Appelhoff`_)
 
-- Add :func:`mne.warp_montage_volume` to use an anatomical volume with channel locations at high or low insensities to apply a symmetric diffeomorphic registration to a template (e.g. fsaverage) space (:gh:`9544` by `Alex Rockhill`_)
+- Add ``mne.warp_montage_volume`` to use an anatomical volume with channel locations at high or low insensities to apply a symmetric diffeomorphic registration to a template (e.g. fsaverage) space (:gh:`9544` by `Alex Rockhill`_)
 
 - :meth:`mne.io.Raw.get_data`, :meth:`mne.Epochs.get_data`, and :meth:`mne.Evoked.get_data` methods now have ``tmin`` and ``tmax`` parameters to specify data window in seconds to get (:gh:`9556` by `Stefan Appelhoff`_)
 

doc/changes/latest.inc

@@ -38,6 +38,7 @@ Enhancements
 - Add a video to :ref:`tut-freesurfer-mne` of a brain inflating from the pial surface to aid in understanding the inflated brain (:gh:`11440` by `Alex Rockhill`_)
 - Add automatic projection of sEEG contact onto the inflated surface for :meth:`mne.viz.Brain.add_sensors` (:gh:`11436` by `Alex Rockhill`_)
 - Allow an image with intracranial electrode contacts (e.g. computed tomography) to be used without the freesurfer recon-all surfaces to locate contacts so that it doesn't have to be downsampled to freesurfer dimensions (for microelectrodes) and show an example :ref:`ex-ieeg-micro` with :func:`mne.transforms.apply_volume_registration_points` added to aid this transform (:gh:`11567` by `Alex Rockhill`_)
+- Use new :meth:`dipy.workflows.align.DiffeomorphicMap.transform_points` to transform a montage of intracranial contacts more efficiently (:gh:`11572` by `Alex Rockhill`_)
 
 Bugs
 ~~~~
@@ -67,3 +68,4 @@ API changes
 ~~~~~~~~~~~
 - Deprecate arguments ``kind`` and ``path`` from :func:`mne.channels.read_layout` in favor of a common argument ``fname`` (:gh:`11500` by `Mathieu Scheltienne`_)
 - Change ``aligned_ct`` positional argument in :func:`mne.gui.locate_ieeg` to ``base_image`` to reflect that this can now be used with unaligned images (:gh:`11567` by `Alex Rockhill`_)
+- ``mne.warp_montage_volume`` was deprecated in favor of :func:`mne.preprocessing.ieeg.warp_montage` (acts directly on points instead of using an intermediate volume) and :func:`mne.preprocessing.ieeg.make_montage_volume` (which makes a volume of ieeg contact locations which can still be useful) (:gh:`11572` by `Alex Rockhill`_)

doc/mri.rst

@@ -35,5 +35,4 @@ Step by step instructions for using :func:`gui.coregistration`:
    transforms.apply_volume_registration_points
    transforms.compute_volume_registration
    vertex_to_mni
-   warp_montage_volume
    coreg.Coregistration

doc/preprocessing.rst

@@ -138,6 +138,8 @@ Projections:
    :toctree: generated/
 
    project_sensors_onto_brain
+   make_montage_volume
+   warp_montage
 
 EEG referencing:
 

mne/preprocessing/ieeg/__init__.py

@@ -5,3 +5,4 @@
 # License: BSD-3-Clause
 
 from ._projection import project_sensors_onto_brain
+from ._volume import make_montage_volume, warp_montage

mne/preprocessing/ieeg/_volume.py

@@ -0,0 +1,219 @@
+# Authors: Alex Rockhill <aprockhill@mailbox.org>
+#
+# License: BSD-3-Clause
+
+import numpy as np
+
+from ...channels import DigMontage, make_dig_montage
+from ...surface import _voxel_neighbors
+from ...transforms import apply_trans, _frame_to_str, Transform
+from ...utils import (verbose, warn, _pl, _validate_type, _require_version,
+                      _check_option)
+
+
+@verbose
+def warp_montage(montage, moving, static, reg_affine, sdr_morph, verbose=None):
+    """Warp a montage to a template with image volumes using SDR.
+
+    .. note:: This is likely only applicable for channels inside the brain
+              (intracranial electrodes).
+
+    Parameters
+    ----------
+    montage : instance of mne.channels.DigMontage
+        The montage object containing the channels.
+    %(moving)s
+    %(static)s
+    %(reg_affine)s
+    %(sdr_morph)s
+    %(verbose)s
+
+    Returns
+    -------
+    montage_warped : mne.channels.DigMontage
+        The modified montage object containing the channels.
+    """
+    _require_version('nibabel', 'warp montage', '2.1.0')
+    _require_version('dipy', 'warping points using SDR', '1.6.0')
+
+    from nibabel import MGHImage
+    from nibabel.spatialimages import SpatialImage
+    from dipy.align.imwarp import DiffeomorphicMap
+
+    _validate_type(moving, SpatialImage, 'moving')
+    _validate_type(static, SpatialImage, 'static')
+    _validate_type(reg_affine, np.ndarray, 'reg_affine')
+    _check_option('reg_affine.shape', reg_affine.shape, ((4, 4),))
+    _validate_type(sdr_morph, (DiffeomorphicMap, None), 'sdr_morph')
+    _validate_type(montage, DigMontage, 'montage')
+
+    moving_mgh = MGHImage(np.array(moving.dataobj).astype(np.float32),
+                          moving.affine)
+    static_mgh = MGHImage(np.array(static.dataobj).astype(np.float32),
+                          static.affine)
+    del moving, static
+
+    # get montage channel coordinates
+    ch_dict = montage.get_positions()
+    if ch_dict['coord_frame'] != 'mri':
+        bad_coord_frames = np.unique([d['coord_frame'] for d in montage.dig])
+        bad_coord_frames = ', '.join([
+            _frame_to_str[cf] if cf in _frame_to_str else str(cf)
+            for cf in bad_coord_frames])
+        raise RuntimeError('Coordinate frame not supported, expected '
+                           f'"mri", got {bad_coord_frames}')
+    ch_names = list(ch_dict['ch_pos'].keys())
+    ch_coords = np.array([ch_dict['ch_pos'][name] for name in ch_names])
+
+    ch_coords = apply_trans(  # convert to moving voxel space
+        np.linalg.inv(moving_mgh.header.get_vox2ras_tkr()), ch_coords * 1000)
+    # next, to moving scanner RAS
+    ch_coords = apply_trans(moving_mgh.header.get_vox2ras(), ch_coords)
+
+    # now, apply reg_affine
+    ch_coords = apply_trans(Transform(  # to static ras
+        fro='ras', to='ras', trans=np.linalg.inv(reg_affine)), ch_coords)
+
+    # now, apply SDR morph
+    if sdr_morph is not None:
+        ch_coords = sdr_morph.transform_points(
+            ch_coords, sdr_morph.domain_grid2world,
+            sdr_morph.domain_world2grid)
+
+    # back to voxels but now for the static image
+    ch_coords = apply_trans(np.linalg.inv(static_mgh.header.get_vox2ras()),
+                            ch_coords)
+
+    # finally, back to surface RAS
+    ch_coords = apply_trans(static_mgh.header.get_vox2ras_tkr(),
+                            ch_coords) / 1000
+
+    # make warped montage
+    montage_warped = make_dig_montage(
+        dict(zip(ch_names, ch_coords)), coord_frame='mri')
+    return montage_warped
+
+
+def _warn_missing_chs(info, dig_image, after_warp=False):
+    """Warn that channels are missing."""
+    # ensure that each electrode contact was marked in at least one voxel
+    missing = set(np.arange(1, len(info.ch_names) + 1)).difference(
+        set(np.unique(np.array(dig_image.dataobj))))
+    missing_ch = [info.ch_names[idx - 1] for idx in missing]
+    if missing_ch:
+        warn(f'Channel{_pl(missing_ch)} '
+             f'{", ".join(repr(ch) for ch in missing_ch)} not assigned '
+             'voxels ' +
+             (f' after applying {after_warp}' if after_warp else ''))
+
+
+@verbose
+def make_montage_volume(montage, base_image, thresh=0.5, max_peak_dist=1,
+                        voxels_max=100, use_min=False, verbose=None):
+    """Make a volume from intracranial electrode contact locations.
+
+    Find areas of the input volume with intensity greater than
+    a threshold surrounding local extrema near the channel location.
+    Monotonicity from the peak is enforced to prevent channels
+    bleeding into each other.
+
+    Parameters
+    ----------
+    montage : instance of mne.channels.DigMontage
+        The montage object containing the channels.
+    base_image : path-like | nibabel.spatialimages.SpatialImage
+        Path to a volumetric scan (e.g. CT) of the subject. Can be in any
+        format readable by nibabel. Can also be a nibabel image object.
+        Local extrema (max or min) should be nearby montage channel locations.
+    thresh : float
+        The threshold relative to the peak to determine the size
+        of the sensors on the volume.
+    max_peak_dist : int
+        The number of voxels away from the channel location to
+        look in the ``image``. This will depend on the accuracy of
+        the channel locations, the default (one voxel in all directions)
+        will work only with localizations that are that accurate.
+    voxels_max : int
+        The maximum number of voxels for each channel.
+    use_min : bool
+        Whether to hypointensities in the volume as channel locations.
+        Default False uses hyperintensities.
+    %(verbose)s
+
+    Returns
+    -------
+    elec_image : nibabel.spatialimages.SpatialImage
+        An image in Freesurfer surface RAS space with voxel values
+        corresponding to the index of the channel. The background
+        is 0s and this index starts at 1.
+    """
+    _require_version('nibabel', 'montage volume', '2.1.0')
+    import nibabel as nib
+
+    _validate_type(montage, DigMontage, 'montage')
+    _validate_type(base_image, nib.spatialimages.SpatialImage, 'base_image')
+    _validate_type(thresh, float, 'thresh')
+    if thresh < 0 or thresh >= 1:
+        raise ValueError(f'`thresh` must be between 0 and 1, got {thresh}')
+    _validate_type(max_peak_dist, int, 'max_peak_dist')
+    _validate_type(voxels_max, int, 'voxels_max')
+    _validate_type(use_min, bool, 'use_min')
+
+    # load image and make sure it's in surface RAS
+    if not isinstance(base_image, nib.spatialimages.SpatialImage):
+        base_image = nib.load(base_image)
+
+    base_image_mgh = nib.MGHImage(
+        np.array(base_image.dataobj).astype(np.float32), base_image.affine)
+    del base_image
+
+    # get montage channel coordinates
+    ch_dict = montage.get_positions()
+    if ch_dict['coord_frame'] != 'mri':
+        bad_coord_frames = np.unique([d['coord_frame'] for d in montage.dig])
+        bad_coord_frames = ', '.join([
+            _frame_to_str[cf] if cf in _frame_to_str else str(cf)
+            for cf in bad_coord_frames])
+        raise RuntimeError('Coordinate frame not supported, expected '
+                           f'"mri", got {bad_coord_frames}')
+
+    ch_names = list(ch_dict['ch_pos'].keys())
+    ch_coords = np.array([ch_dict['ch_pos'][name] for name in ch_names])
+
+    # convert to voxel space
+    ch_coords = apply_trans(
+        np.linalg.inv(base_image_mgh.header.get_vox2ras_tkr()),
+        ch_coords * 1000)
+
+    # take channel coordinates and use the image to transform them
+    # into a volume where all the voxels over a threshold nearby
+    # are labeled with an index
+    image_data = np.array(base_image_mgh.dataobj)
+    if use_min:
+        image_data *= -1
+    elec_image = np.zeros(base_image_mgh.shape, dtype=int)
+    for i, ch_coord in enumerate(ch_coords):
+        if np.isnan(ch_coord).any():
+            continue
+        # this looks up to a voxel away, it may be marked imperfectly
+        volume = _voxel_neighbors(ch_coord, image_data, thresh=thresh,
+                                  max_peak_dist=max_peak_dist,
+                                  voxels_max=voxels_max)
+        for voxel in volume:
+            if elec_image[voxel] != 0:
+                # some voxels ambiguous because the contacts are bridged on
+                # the image so assign the voxel to the nearest contact location
+                dist_old = np.sqrt(
+                    (ch_coords[elec_image[voxel] - 1] - voxel)**2).sum()
+                dist_new = np.sqrt((ch_coord - voxel)**2).sum()
+                if dist_new < dist_old:
+                    elec_image[voxel] = i + 1
+            else:
+                elec_image[voxel] = i + 1
+
+    # assemble the volume
+    elec_image = nib.spatialimages.SpatialImage(
+        elec_image, base_image_mgh.affine)
+    _warn_missing_chs(montage, elec_image, after_warp=False)
+
+    return elec_image

mne/preprocessing/ieeg/tests/test_projection.py

@@ -1,8 +1,5 @@
 # -*- coding: utf-8 -*-
-"""Test the compute_current_source_density function.
-
-For each supported file format, implement a test.
-"""
+"""Test the ieeg projection functions."""
 # Authors: Alex Rockhill <aprockhill@mailbox.org>
 #
 # License: BSD-3-Clause

mne/preprocessing/ieeg/tests/test_volume.py

@@ -0,0 +1,105 @@
+# -*- coding: utf-8 -*-
+"""Test ieeg volume functions."""
+# Authors: Alex Rockhill <aprockhill@mailbox.org>
+#
+# License: BSD-3-Clause
+
+import numpy as np
+import pytest
+
+from mne.coreg import get_mni_fiducials
+from mne.channels import make_dig_montage
+from mne.datasets import testing
+from mne.preprocessing.ieeg import make_montage_volume, warp_montage
+from mne.transforms import apply_trans, compute_volume_registration
+
+data_path = testing.data_path(download=False)
+subjects_dir = data_path / "subjects"
+
+
+@pytest.mark.slowtest
+@testing.requires_testing_data
+def test_warp_montage():
+    """Test warping an montage based on intracranial electrode positions."""
+    nib = pytest.importorskip('nibabel')
+    pytest.importorskip('dipy')
+    subject_brain = nib.load(subjects_dir / "sample" / "mri" / "brain.mgz")
+    template_brain = nib.load(subjects_dir / "fsaverage" / "mri" / "brain.mgz")
+    zooms = dict(translation=10, rigid=10, sdr=10)
+    reg_affine, sdr_morph = compute_volume_registration(
+        subject_brain, template_brain, zooms=zooms,
+        niter=[3, 3, 3],
+        pipeline=('translation', 'rigid', 'sdr'))
+    # make an info object with three channels with positions
+    ch_coords = np.array([[-8.7040273, 17.99938754, 10.29604017],
+                          [-14.03007764, 19.69978401, 12.07236939],
+                          [-21.1130506, 21.98310911, 13.25658887]])
+    ch_pos = dict(zip(['1', '2', '3'], ch_coords / 1000))  # mm -> m
+    lpa, nasion, rpa = get_mni_fiducials('sample', subjects_dir)
+    montage = make_dig_montage(ch_pos, lpa=lpa['r'], nasion=nasion['r'],
+                               rpa=rpa['r'], coord_frame='mri')
+    montage_warped = warp_montage(
+        montage, subject_brain, template_brain, reg_affine, sdr_morph)
+    # checked with nilearn plot from `tut-ieeg-localize`
+    # check montage in surface RAS
+    ground_truth_warped = np.array([[-0.009, -0.00133333, -0.033],
+                                    [-0.01445455, 0.00127273, -0.03163636],
+                                    [-0.022, 0.00285714, -0.031]])
+    for i, d in enumerate(montage_warped.dig):
+        assert np.linalg.norm(  # off by less than 1 cm
+            d['r'] - ground_truth_warped[i]) < 0.01
+
+    bad_montage = montage.copy()
+    for d in bad_montage.dig:
+        d['coord_frame'] = 99
+    with pytest.raises(RuntimeError, match='Coordinate frame not supported'):
+        warp_montage(bad_montage, subject_brain, template_brain, reg_affine,
+                     sdr_morph)
+
+
+@pytest.mark.slowtest
+@testing.requires_testing_data
+def test_make_montage_volume():
+    """Test making a montage image based on intracranial electrodes."""
+    nib = pytest.importorskip('nibabel')
+    pytest.importorskip('dipy')
+    subject_brain = nib.load(subjects_dir / "sample" / "mri" / "brain.mgz")
+    # make an info object with three channels with positions
+    ch_coords = np.array([[-8.7040273, 17.99938754, 10.29604017],
+                          [-14.03007764, 19.69978401, 12.07236939],
+                          [-21.1130506, 21.98310911, 13.25658887]])
+    ch_pos = dict(zip(['1', '2', '3'], ch_coords / 1000))  # mm -> m
+    lpa, nasion, rpa = get_mni_fiducials('sample', subjects_dir)
+    montage = make_dig_montage(ch_pos, lpa=lpa['r'], nasion=nasion['r'],
+                               rpa=rpa['r'], coord_frame='mri')
+    # make fake image based on the info
+    CT_data = np.zeros(subject_brain.shape)
+    # convert to voxels
+    ch_coords_vox = apply_trans(
+        np.linalg.inv(subject_brain.header.get_vox2ras_tkr()), ch_coords)
+    for (x, y, z) in ch_coords_vox.round().astype(int):
+        # make electrode contact hyperintensities
+        # first, make the surrounding voxels high intensity
+        CT_data[x - 1:x + 2, y - 1:y + 2, z - 1:z + 2] = 500
+        # then, make the center even higher intensity
+        CT_data[x, y, z] = 1000
+    CT = nib.Nifti1Image(CT_data, subject_brain.affine)
+
+    elec_image = make_montage_volume(montage, CT, thresh=0.25)
+    elec_image_data = np.array(elec_image.dataobj)
+
+    # check elec image, center should be no more than half a voxel away
+    for i in range(len(montage.ch_names)):
+        assert np.linalg.norm(
+            np.array(np.where(elec_image_data == i + 1)
+                     ).mean(axis=1) - ch_coords_vox[i]) < 0.5
+
+    # test inputs
+    with pytest.raises(ValueError, match='`thresh` must be between 0 and 1'):
+        make_montage_volume(montage, CT, thresh=11.)
+
+    bad_montage = montage.copy()
+    for d in bad_montage.dig:
+        d['coord_frame'] = 99
+    with pytest.raises(RuntimeError, match='Coordinate frame not supported'):
+        make_montage_volume(bad_montage, CT)