[DOC, MRG] Fix CT alignment instructions (mne-tools#10298)

* fix CT alignment instructions * add import * fix narration in code block * fix narrative for computationally inefficient version * Eric version * add doc reference * Eric comment clarification * Eric comment * FIX: nilearn * fix conf * fix conf try again * try another fix Co-authored-by: Eric Larson <[email protected]>
jevenzh · Feb 9, 2022 · c3e812b · c3e812b
1 parent e6b5367
commit c3e812b
Show file tree

Hide file tree

Showing 6 changed files with 103 additions and 23 deletions.
diff --git a/doc/conf.py b/doc/conf.py
@@ -832,6 +832,7 @@ def reset_warnings(gallery_conf, fname):
                 "default value of type 'dict' in an Any trait will",  # traits
                 'rcParams is deprecated',  # PyVista rcParams -> global_theme
                 'to mean no clipping',
+                r'the `scipy\.ndimage.*` namespace is deprecated',  # Dipy
                 '`np.MachAr` is deprecated',  # Numba
                 ):
         warnings.filterwarnings(  # deal with other modules having bad imports

diff --git a/doc/mri.rst b/doc/mri.rst
@@ -25,6 +25,7 @@ Step by step instructions for using :func:`gui.coregistration`:
    head_to_mni
    head_to_mri
    read_freesurfer_lut
+   read_lta
    read_talxfm
    scale_mri
    scale_bem

diff --git a/mne/__init__.py b/mne/__init__.py
@@ -46,7 +46,7 @@
                     find_stim_steps, AcqParserFIF)
 from ._freesurfer import (head_to_mni, head_to_mri, read_talxfm,
                           get_volume_labels_from_aseg, read_freesurfer_lut,
-                          vertex_to_mni)
+                          vertex_to_mni, read_lta)
 from .forward import (read_forward_solution, apply_forward, apply_forward_raw,
                       average_forward_solutions, Forward,
                       write_forward_solution, make_forward_solution,

diff --git a/mne/_freesurfer.py b/mne/_freesurfer.py
@@ -463,6 +463,28 @@ def _ensure_image_in_surface_RAS(image, subject, subjects_dir):
     return image  # returns MGH image for header
 
 
+@verbose
+def read_lta(fname, verbose=None):
+    """Read a Freesurfer linear transform array file.
+
+    Parameters
+    ----------
+    fname : str | None
+        The transform filename.
+    %(verbose)s
+
+    Returns
+    -------
+    affine : ndarray
+        The affine transformation described by the lta file.
+    """
+    _validate_type(fname, ('path-like', None), 'fname')
+    _check_fname(fname, 'read', must_exist=True)
+    with open(fname, 'r') as fid:
+        affine = np.loadtxt(fid.readlines()[5:9])
+    return affine
+
+
 @verbose
 def read_talxfm(subject, subjects_dir=None, verbose=None):
     """Compute MRI-to-MNI transform from FreeSurfer talairach.xfm file.

diff --git a/mne/tests/test_freesurfer.py b/mne/tests/test_freesurfer.py
@@ -9,7 +9,7 @@
                  read_talxfm, read_freesurfer_lut,
                  get_volume_labels_from_aseg)
 from mne.datasets import testing
-from mne._freesurfer import _get_mgz_header, _check_subject_dir
+from mne._freesurfer import _get_mgz_header, _check_subject_dir, read_lta
 from mne.transforms import apply_trans, _get_trans
 from mne.utils import requires_nibabel
 
@@ -96,6 +96,44 @@ def test_vertex_to_mni_fs_nibabel(monkeypatch):
     assert_allclose(coords, coords_2, atol=0.1)
 
 
+def test_read_lta(tmp_path):
+    """Test reading a Freesurfer linear transform array file."""
+    with open(op.join(tmp_path, 'test.lta'), 'w') as fid:
+        fid.write("""type      = 0 # LINEAR_VOX_TO_VOX
+                     nxforms   = 1
+                     mean      = 0.0000 0.0000 0.0000
+                     sigma     = 1.0000
+                     1 4 4
+                     0.99221027 -0.05494503  0.11180324 -3.84350586
+                     0.05233596  0.99828744 0.02614108 -9.77523804
+                     -0.11304809 -0.02008611 0.99338663 15.25457001
+                     0 0 0 1
+                     src volume info
+                     valid = 1  # volume info valid
+                     filename = tmp.mgz
+                     volume = 256 256 256
+                     voxelsize = 1 1 1
+                     xras   = -1 0 0
+                     yras   = 0 0 -1
+                     zras   = 0 1 0
+                     cras   = -1.19374 -3.31686 3.25835
+                     dst volume info
+                     valid = 1  # volume info valid
+                     filename = tmp.mgz
+                     volume = 256 256 256
+                     voxelsize = 1 1 1
+                     xras   = -1 0 0
+                     yras   = 0 0 -1
+                     zras   = 0 1 0
+                     cras   = -1.19374 -3.31686 3.25835)""")
+    assert_array_equal(
+        read_lta(op.join(tmp_path, 'test.lta')),
+        np.array([[0.99221027, -0.05494503, 0.11180324, -3.84350586],
+                  [0.05233596, 0.99828744, 0.02614108, -9.77523804],
+                  [-0.11304809, -0.02008611, 0.99338663, 15.25457001],
+                  [0., 0., 0., 1.]]))
+
+
 @testing.requires_testing_data
 @requires_nibabel()
 @pytest.mark.parametrize('fname', [

diff --git a/tutorials/clinical/10_ieeg_localize.py b/tutorials/clinical/10_ieeg_localize.py
@@ -187,7 +187,7 @@ def plot_overlay(image, compare, title, thresh=None):
 # here::
 #
 #    reg_affine, _ = mne.transforms.compute_volume_registration(
-#         CT_orig, T1, pipeline='rigids')
+#         CT_orig, T1, pipeline='rigids', zooms=dict(translation=5)))
 #
 # And instead we just hard-code the resulting 4x4 matrix:
 
@@ -202,26 +202,44 @@ def plot_overlay(image, compare, title, thresh=None):
 
 # %%
 # .. note::
-#     Alignment failures sometimes occur which requires manual alignment.
-#     This can be done using Freesurfer's ``freeview`` to align manually
-#
-#         - Load the two scans from the command line using
-#           ``freeview $MISC_PATH/seeg/sample_seeg/mri/T1.mgz
-#           $MISC_PATH/seeg/sample_seeg_CT.mgz``
-#         - Navigate to the upper toolbar, go to
-#           :menuselection:`Tools --> Transform Volume`
-#         - Use the rotation and translation slide bars to align the CT
-#           to the MR (be sure to have the CT selected in the upper left menu)
-#         - Save the modified volume using the ``Save Volume As...`` button
-#         - Resample to the T1 shape and affine using::
-#
-#               CT_aligned_pre = nib.load(op.join(misc_path, 'seeg',
-#                                                 'sample_seeg_CT_aligned.mgz'))
-#               CT_aligned = resample(
-#                   moving=np.asarray(CT_aligned_pre.dataobj),
-#                   static=np.asarray(T1.dataobj),
-#                   moving_affine=CT_aligned_pre.affine,
-#                   static_affine=T1.affine)
+#     Alignment failures sometimes occur which requires manual pre-alignment.
+#     Freesurfer's ``freeview`` can be used to to align manually
+#
+#     .. code-block:: bash
+#
+#         $ freeview $MISC_PATH/seeg/sample_seeg/mri/T1.mgz \
+#            $MISC_PATH/seeg/sample_seeg_CT.mgz:colormap=heat:opacity=0.6
+#
+#     - Navigate to the upper toolbar, go to
+#       :menuselection:`Tools --> Transform Volume`
+#     - Use the rotation and translation slide bars to align the CT
+#       to the MR (be sure to have the CT selected in the upper left menu)
+#     - Save the linear transform array (lta) file using the ``Save Reg...``
+#       button
+#
+#     Since we really require as much precision as possible for the
+#     alignment, we should rerun the algorithm starting with the manual
+#     alignment. This time, we just want to skip to the most exact rigid
+#     alignment, without smoothing, since the manual alignment is already
+#     very close.
+#
+#     .. code-block:: python
+#
+#         from dipy.align import affine_registration
+#         # load transform
+#         manual_reg_affine_vox = mne.read_lta(op.join(  # the path used above
+#             misc_path, 'seeg', 'sample_seeg_CT_aligned_manual.mgz.lta'))
+#         # convert from vox->vox to ras->ras
+#         manual_reg_affine = \
+#             CT_orig.affine @ np.linalg.inv(manual_reg_affine_vox) \
+#             @ np.linalg.inv(T1.affine)
+#         CT_aligned_fix_img = affine_registration(
+#             moving=np.array(CT_orig.dataobj), static=np.array(T1.dataobj),
+#             moving_affine=CT_orig.affine, static_affine=T1.affine,
+#             pipeline=['rigid'], starting_affine=manual_reg_affine,
+#             level_iters=[100], sigmas=[0], factors=[1])[0]
+#         CT_aligned = nib.MGHImage(
+#             CT_aligned_fix_img.astype(np.float32), T1.affine)
 #
 #     The rest of the tutorial can then be completed using ``CT_aligned``
 #     from this point on.