diff --git a/doc/conf.py b/doc/conf.py index 264945f039f..8c1ab65ddff 100644 --- 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 index 9d4fa797c4e..483c14277ce 100644 --- 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 index b5624a68628..82a9fa3926a 100644 --- 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 index 13ff80eecf3..7b5e33f808d 100644 --- 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 index c2583c36812..4f154b24cbf 100644 --- 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 index e25d3482bf5..9d830656876 100644 --- 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.