MRG, ENH: Abstracted volume registration again

doc/changes/latest.inc

@@ -73,7 +73,9 @@ Enhancements
 
 - Show multiple colors and linestyles for excluded components with :class:`mne.Evoked` in :meth:`mne.preprocessing.ICA.plot_sources` (:gh:`9444` by `Martin Schulz`_)
 
-- Added tutorial for how to processes image (CT and MR) files in order to localize electrode contacts for intracranial recordings :ref:`tut-ieeg-localize` (:gh`9484` by `Alex Rockhill`_)
+- Add functions for aligning MRI and CT data `mne.transforms.compute_volume_registration` and `mne.transforms.apply_volume_registration` (:gh:`9503` by `Alex Rockhill`_ and `Eric Larson`_)
+
+- Add tutorial for how to processes image (CT and MR) files in order to localize electrode contacts for intracranial recordings :ref:`tut-ieeg-localize` (:gh`9484` by `Alex Rockhill`_)
 
 - Add support for colormap normalization in :func:`mne.viz.plot_topomap` (:gh:`9468` by `Clemens Brunner`_)
 
@@ -95,8 +97,6 @@ Bugs
 
 - Fix bug with `mne.io.read_raw_fieldtrip` and `mne.read_epochs_fieldtrip` where channel positions were not set properly (:gh:`9447` by `Eric Larson`_)
 
-- Fix bug with `mne.Epochs.plot` where event lines were misplaced if epochs were overlapping in time (:gh:`9505` by `Daniel McCloy`_)
-
 - :func:`mne.concatenate_raws` now raises an exception if ``raw.info['dev_head_t']`` differs between files. This behavior can be controlled using the new ``on_mismatch`` parameter (:gh:`9438` by `Richard Höchenberger`_)
 
 - Fixed bug in :meth:`mne.Epochs.drop_bad` where subsequent rejections failed if they only specified thresholds for a subset of the channel types used in a previous rejection (:gh:`9485` by `Richard Höchenberger`_).
@@ -111,4 +111,4 @@ Bugs
 
 API changes
 ~~~~~~~~~~~
-- Nothing yet
+- In `mne.compute_source_morph`, the ``niter_affine`` and ``niter_sdr`` parameters have been replaced by ``niter`` and ``pipeline`` parameters for more consistent and finer-grained control of registration/warping steps and iteration (:gh:`9505` by `Alex Rockhill`_ and `Eric Larson`_)

doc/mri.rst

@@ -18,9 +18,11 @@ Step by step instructions for using :func:`gui.coregistration`:
    gui.coregistration
    gui.fiducials
    create_default_subject
-   marching_cubes
    scale_mri
    scale_bem
    scale_labels
    scale_source_space
+   surface.marching_cubes
+   transforms.apply_volume_registration
+   transforms.compute_volume_registration
    voxel_neighbors

mne/__init__.py

@@ -63,7 +63,7 @@
                               extract_label_time_course, stc_near_sensors)
 from .surface import (read_surface, write_surface, decimate_surface, read_tri,
                       get_head_surf, get_meg_helmet_surf, dig_mri_distances,
-                      marching_cubes, voxel_neighbors)
+                      voxel_neighbors)
 from .morph_map import read_morph_map
 from .morph import (SourceMorph, read_source_morph, grade_to_vertices,
                     compute_source_morph)

mne/datasets/utils.py

@@ -254,7 +254,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
     path = _get_path(path, key, name)
     # To update the testing or misc dataset, push commits, then make a new
     # release on GitHub. Then update the "releases" variable:
-    releases = dict(testing='0.121', misc='0.11')
+    releases = dict(testing='0.121', misc='0.13')
     # And also update the "md5_hashes['testing']" variable below.
     # To update any other dataset, update the data archive itself (upload
     # an updated version) and update the md5 hash.
@@ -345,7 +345,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
             bst_raw='fa2efaaec3f3d462b319bc24898f440c',
             bst_resting='70fc7bf9c3b97c4f2eab6260ee4a0430'),
         fake='3194e9f7b46039bb050a74f3e1ae9908',
-        misc='4c07fe4e3a068301417362e835c9af59',
+        misc='ab4be0b325f5896a28165ad177bda338',
         sample='12b75d1cb7df9dfb4ad73ed82f61094f',
         somato='32fd2f6c8c7eb0784a1de6435273c48b',
         spm='9f43f67150e3b694b523a21eb929ea75',

mne/defaults.py

@@ -99,7 +99,14 @@
         alpha=None, resolution=1., surface_alpha=None, blending='mip',
         silhouette_alpha=None, silhouette_linewidth=2.),
     prefixes={'': 1e0, 'd': 1e1, 'c': 1e2, 'm': 1e3, 'µ': 1e6, 'u': 1e6,
-              'n': 1e9, 'p': 1e12, 'f': 1e15}
+              'n': 1e9, 'p': 1e12, 'f': 1e15},
+    transform_zooms=dict(
+        translation=None, rigid=None, affine=None, sdr=None),
+    transform_niter=dict(
+        translation=(100, 100, 10),
+        rigid=(100, 100, 10),
+        affine=(100, 100, 10),
+        sdr=(5, 5, 3)),
 )
 
 

mne/morph.py

@@ -17,10 +17,9 @@
     _get_ico_tris)
 from .source_space import SourceSpaces, _ensure_src, _grid_interp
 from .surface import mesh_edges, read_surface, _compute_nearest
-from .transforms import _angle_between_quats, rot_to_quat
 from .utils import (logger, verbose, check_version, get_subjects_dir,
                     warn as warn_, fill_doc, _check_option, _validate_type,
-                    BunchConst, wrapped_stdout, _check_fname, warn,
+                    BunchConst, _check_fname, warn,
                     _ensure_int, ProgressBar, use_log_level)
 from .externals.h5io import read_hdf5, write_hdf5
 
@@ -1006,100 +1005,20 @@ def _interpolate_data(stc, morph, mri_resolution, mri_space, output):
 ###############################################################################
 # Morph for VolSourceEstimate
 
-def _compute_r2(a, b):
-    return 100 * (a.ravel() @ b.ravel()) / \
-        (np.linalg.norm(a) * np.linalg.norm(b))
-
-
 def _compute_morph_sdr(mri_from, mri_to, niter_affine, niter_sdr, zooms):
     """Get a matrix that morphs data from one subject to another."""
-    with np.testing.suppress_warnings():
-        from dipy.align import imaffine, imwarp, metrics, transforms
-    from dipy.align.reslice import reslice
-
-    logger.info('Computing nonlinear Symmetric Diffeomorphic Registration...')
-
-    # reslice mri_from to zooms
-    mri_from_orig = mri_from
-    mri_from, mri_from_affine = reslice(
-        _get_img_fdata(mri_from_orig), mri_from_orig.affine,
-        mri_from_orig.header.get_zooms()[:3], zooms)
-
-    # reslice mri_to to zooms
-    mri_to, affine = reslice(
-        _get_img_fdata(mri_to), mri_to.affine,
-        mri_to.header.get_zooms()[:3], zooms)
-
-    mri_to /= mri_to.max()
-    mri_from /= mri_from.max()  # normalize
-
-    # compute center of mass
-    c_of_mass = imaffine.transform_centers_of_mass(
-        mri_to, affine, mri_from, mri_from_affine)
-
-    # set up Affine Registration
-    affreg = imaffine.AffineRegistration(
-        metric=imaffine.MutualInformationMetric(nbins=32),
-        level_iters=list(niter_affine),
-        sigmas=[3.0, 1.0, 0.0],
-        factors=[4, 2, 1])
-
-    # translation
-    logger.info('Optimizing translation:')
-    with wrapped_stdout(indent='    ', cull_newlines=True):
-        translation = affreg.optimize(
-            mri_to, mri_from, transforms.TranslationTransform3D(), None,
-            affine, mri_from_affine, starting_affine=c_of_mass.affine)
-    mri_from_to = translation.transform(mri_from)
-    dist = np.linalg.norm(translation.affine[:3, 3])
-    logger.info(f'    Translation: {dist:6.1f} mm')
-    logger.info(f'    R²:          {_compute_r2(mri_to, mri_from_to):6.1f}%')
-
-    # rigid body transform (translation + rotation)
-    logger.info('Optimizing rigid-body:')
-    with wrapped_stdout(indent='    ', cull_newlines=True):
-        rigid = affreg.optimize(
-            mri_to, mri_from, transforms.RigidTransform3D(), None,
-            affine, mri_from_affine, starting_affine=translation.affine)
-    mri_from_to = rigid.transform(mri_from)
-    dist = np.linalg.norm(rigid.affine[:3, 3])
-    angle = np.rad2deg(_angle_between_quats(
-        np.zeros(3), rot_to_quat(rigid.affine[:3, :3])))
-
-    logger.info(f'    Translation: {dist:6.1f} mm')
-    logger.info(f'    Rotation:    {angle:6.1f}°')
-    logger.info(f'    R²:          {_compute_r2(mri_to, mri_from_to):6.1f}%')
-
-    # affine transform (translation + rotation + scaling)
-    logger.info('Optimizing full affine:')
-    with wrapped_stdout(indent='    ', cull_newlines=True):
-        pre_affine = affreg.optimize(
-            mri_to, mri_from, transforms.AffineTransform3D(), None,
-            affine, mri_from_affine, starting_affine=rigid.affine)
-    mri_from_to = pre_affine.transform(mri_from)
-    logger.info(
-        f'    R²:          {_compute_r2(mri_to, mri_from_to):6.1f}%')
-
-    # SDR
-    shape = tuple(pre_affine.domain_shape)
-    if len(niter_sdr):
-        sdr = imwarp.SymmetricDiffeomorphicRegistration(
-            metrics.CCMetric(3), list(niter_sdr))
-        logger.info('Optimizing SDR:')
-        with wrapped_stdout(indent='    ', cull_newlines=True):
-            sdr_morph = sdr.optimize(mri_to, pre_affine.transform(mri_from))
-        assert shape == tuple(sdr_morph.domain_shape)  # should be tuple of int
-        mri_from_to = sdr_morph.transform(mri_from_to)
-        logger.info(
-            f'    R²:          {_compute_r2(mri_to, mri_from_to):6.1f}%')
-    else:
-        sdr_morph = None
-
-    _debug_img(mri_from_orig.dataobj, mri_from_orig.affine, 'From')
-    _debug_img(mri_from, affine, 'From-reslice')
-    _debug_img(mri_from_to, affine, 'From-reslice')
-    _debug_img(mri_to, affine, 'To-reslice')
-    return shape, zooms, affine, pre_affine, sdr_morph
+    from .transforms import _compute_volume_registration
+    from dipy.align.imaffine import AffineMap
+    pipeline = 'all' if niter_sdr else 'affines'
+    niter = dict(translation=niter_affine, rigid=niter_affine,
+                 affine=niter_affine,
+                 sdr=niter_sdr if niter_sdr else (1,))
+    pre_affine, sdr_morph, to_shape, to_affine, from_shape, from_affine = \
+        _compute_volume_registration(
+            mri_from, mri_to, zooms=zooms, niter=niter, pipeline=pipeline)
+    pre_affine = AffineMap(
+        pre_affine, to_shape, to_affine, from_shape, from_affine)
+    return to_shape, zooms, to_affine, pre_affine, sdr_morph
 
 
 def _compute_morph_matrix(subject_from, subject_to, vertices_from, vertices_to,

mne/tests/test_morph.py

@@ -467,11 +467,11 @@ def test_volume_source_morph_basic(tmpdir):
 @testing.requires_testing_data
 @pytest.mark.parametrize(
     'subject_from, subject_to, lower, upper, dtype, morph_mat', [
-        ('sample', 'fsaverage', 5.9, 6.1, float, False),
+        ('sample', 'fsaverage', 10.0, 10.4, float, False),
         ('fsaverage', 'fsaverage', 0., 0.1, float, False),
         ('sample', 'sample', 0., 0.1, complex, False),
         ('sample', 'sample', 0., 0.1, float, True),  # morph_mat
-        ('sample', 'fsaverage', 10, 12, float, True),  # morph_mat
+        ('sample', 'fsaverage', 7.0, 7.4, float, True),  # morph_mat
     ])
 def test_volume_source_morph_round_trip(
         tmpdir, subject_from, subject_to, lower, upper, dtype, morph_mat,
@@ -554,7 +554,7 @@ def test_volume_source_morph_round_trip(
     # check that power is more or less preserved (labelizing messes with this)
     if morph_mat:
         if subject_to == 'fsaverage':
-            limits = (18, 18.5)
+            limits = (14.0, 14.2)
         else:
             limits = (7, 7.5)
     else:

mne/tests/test_source_estimate.py

@@ -1803,7 +1803,7 @@ def test_scale_morph_labels(kind, scale, monkeypatch, tmpdir):
                 min_, max_ = 0.72, 0.75
             else:
                 # min_, max_ = 0.84, 0.855  # zooms='auto' values
-                min_, max_ = 0.61, 0.62
+                min_, max_ = 0.61, 0.63
             assert min_ < corr <= max_, scale
         else:
             assert_allclose(

mne/tests/test_transforms.py

@@ -8,6 +8,7 @@
 
 import mne
 from mne.datasets import testing
+from mne.fixes import _get_img_fdata
 from mne import read_trans, write_trans
 from mne.io import read_info
 from mne.transforms import (invert_transform, _get_trans,
@@ -21,12 +22,15 @@
                             rotation3d_align_z_axis, _read_fs_xfm,
                             _write_fs_xfm, _quat_real, _fit_matched_points,
                             _quat_to_euler, _euler_to_quat,
-                            _quat_to_affine)
+                            _quat_to_affine, _compute_r2)
+from mne.utils import requires_nibabel, requires_dipy
 
 data_path = testing.data_path(download=False)
 fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-trans.fif')
 fname_eve = op.join(data_path, 'MEG', 'sample',
                     'sample_audvis_trunc_raw-eve.fif')
+subjects_dir = op.join(data_path, 'subjects')
+fname_t1 = op.join(subjects_dir, 'fsaverage', 'mri', 'T1.mgz')
 
 base_dir = op.join(op.dirname(__file__), '..', 'io', 'tests', 'data')
 fname_trans = op.join(base_dir, 'sample-audvis-raw-trans.txt')
@@ -507,3 +511,29 @@ def test_euler(quats):
     quat_rot = quat_to_rot(quats)
     euler_rot = np.array([rotation(*e)[:3, :3] for e in euler])
     assert_allclose(quat_rot, euler_rot, atol=1e-14)
+
+
+@requires_nibabel()
+@requires_dipy()
+@pytest.mark.slowtest
+@testing.requires_testing_data
+def test_volume_registration():
+    """Test volume registration."""
+    import nibabel as nib
+    from dipy.align import resample
+    T1 = nib.load(fname_t1)
+    affine = np.eye(4)
+    affine[0, 3] = 10
+    T1_resampled = resample(moving=T1.get_fdata(),
+                            static=T1.get_fdata(),
+                            moving_affine=T1.affine,
+                            static_affine=T1.affine,
+                            between_affine=np.linalg.inv(affine))
+    for pipeline in ('rigids', ('translation', 'sdr')):
+        reg_affine, sdr_morph = mne.transforms.compute_volume_registration(
+            T1_resampled, T1, pipeline=pipeline, zooms=10, niter=[5])
+        assert_allclose(affine, reg_affine, atol=0.25)
+        T1_aligned = mne.transforms.apply_volume_registration(
+            T1_resampled, T1, reg_affine, sdr_morph)
+        r2 = _compute_r2(_get_img_fdata(T1_aligned), _get_img_fdata(T1))
+        assert 99.9 < r2