Skip to content

Commit cb132a6

Browse files
alexrockhillalexrockhill
committed
working version, lot's of reverting
1 parent a83cdb3 commit cb132a6

7 files changed

Lines changed: 275 additions & 311 deletions

File tree

‎mne/defaults.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -101,7 +101,7 @@
101101
prefixes={'': 1e0, 'd': 1e1, 'c': 1e2, 'm': 1e3, 'µ': 1e6, 'u': 1e6,
102102
'n': 1e9, 'p': 1e12, 'f': 1e15},
103103
transform_zooms=dict(
104-
translation=None, rigid=None, affine=None, sdr=None),
104+
translation=5, rigid=None, affine=None, sdr=None),
105105
transform_niter=dict(
106106
translation=(100, 100, 10),
107107
rigid=(100, 100, 10),

‎mne/morph.py‎

Lines changed: 29 additions & 207 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,6 @@
99
import copy
1010
import numpy as np
1111

12-
from .defaults import _handle_default
1312
from .fixes import _get_img_fdata
1413
from .morph_map import read_morph_map
1514
from .parallel import parallel_func
@@ -18,11 +17,12 @@
1817
_get_ico_tris)
1918
from .source_space import SourceSpaces, _ensure_src, _grid_interp
2019
from .surface import mesh_edges, read_surface, _compute_nearest
21-
from .transforms import _angle_between_quats, rot_to_quat
20+
from .transforms import (compute_volume_registration, _reslice_normalize,
21+
_validate_niter)
2222
from .utils import (logger, verbose, check_version, get_subjects_dir,
2323
warn as warn_, fill_doc, _check_option, _validate_type,
2424
BunchConst, wrapped_stdout, _check_fname, warn,
25-
_ensure_int, ProgressBar, use_log_level)
25+
_ensure_int, ProgressBar, use_log_level, _check_dep)
2626
from .externals.h5io import read_hdf5, write_hdf5
2727

2828

@@ -181,7 +181,7 @@ def compute_source_morph(src, subject_from=None, subject_to='fsaverage',
181181
'sparse morph.')
182182

183183
subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
184-
shape = affine = pre_affine = sdr_morph = morph_mat = None
184+
shape = pre_affine = sdr_morph = morph_mat = None
185185
vertices_to_surf, vertices_to_vol = list(), list()
186186

187187
# niter_affine and niter_sdr deprecations (0.24)
@@ -212,6 +212,8 @@ def compute_source_morph(src, subject_from=None, subject_to='fsaverage',
212212
if kind in ('volume', 'mixed'):
213213
_check_dep(nibabel='2.1.0', dipy='0.10.1')
214214
import nibabel as nib
215+
with np.testing.suppress_warnings():
216+
from dipy.align import imaffine
215217

216218
logger.info('Volume source space(s) present...')
217219

@@ -252,9 +254,13 @@ def compute_source_morph(src, subject_from=None, subject_to='fsaverage',
252254
zooms_src_to = tuple(zooms_src_to)
253255

254256
# pre-compute non-linear morph
255-
zooms = _check_zooms(mri_from, zooms, zooms_src_to)
256-
affine, pre_affine, sdr_morph = _compute_morph_sdr(
257-
mri_from, mri_to, pipeline, niter, zooms)
257+
pre_affine, sdr_morph = compute_volume_registration(
258+
mri_from, mri_to, pipeline, niter,
259+
5. if zooms == 'auto' else zooms)
260+
# convert to affine map
261+
pre_affine = imaffine.AffineMap(
262+
pre_affine, mri_to.shape, mri_to.affine,
263+
mri_from.shape, mri_to.affine)
258264
shape = tuple(pre_affine.domain_shape)
259265

260266
if kind in ('surface', 'mixed'):
@@ -289,7 +295,7 @@ def compute_source_morph(src, subject_from=None, subject_to='fsaverage',
289295
assert len(vertices_to) == len(src_to)
290296
morph = SourceMorph(subject_from, subject_to, kind, zooms,
291297
niter['affine'], niter['sdr'], spacing, smooth, xhemi,
292-
morph_mat, vertices_to, shape, affine,
298+
morph_mat, vertices_to, shape, mri_to.affine,
293299
pre_affine, sdr_morph, src_data, None)
294300
if precompute:
295301
morph.compute_vol_morph_mat()
@@ -700,22 +706,6 @@ def _debug_img(data, affine, title, shape=None):
700706
return
701707

702708

703-
def _check_zooms(mri_from, zooms, zooms_src_to):
704-
# use voxel size of mri_from
705-
if isinstance(zooms, str) and zooms == 'auto':
706-
zooms = zooms_src_to if zooms_src_to is not None else 5.
707-
if zooms is None:
708-
zooms = mri_from.header.get_zooms()[:3]
709-
zooms = np.atleast_1d(zooms).astype(float)
710-
if zooms.shape == (1,):
711-
zooms = np.repeat(zooms, 3)
712-
if zooms.shape != (3,):
713-
raise ValueError('zooms must be None, a singleton, or have shape (3,),'
714-
' got shape %s' % (zooms.shape,))
715-
zooms = tuple(zooms)
716-
return zooms
717-
718-
719709
def _resample_from_to(img, affine, to_vox_map):
720710
# Wrap to dipy for speed, equivalent to:
721711
# from nibabel.processing import resample_from_to
@@ -789,24 +779,13 @@ def read_source_morph(fname):
789779

790780
###############################################################################
791781
# Helper functions for SourceMorph methods
792-
def _check_dep(nibabel='2.1.0', dipy='0.10.1'):
793-
"""Check dependencies."""
794-
for lib, ver in zip(['nibabel', 'dipy'],
795-
[nibabel, dipy]):
796-
passed = True if not ver else check_version(lib, ver)
797-
798-
if not passed:
799-
raise ImportError('%s %s or higher must be correctly '
800-
'installed and accessible from Python' % (lib,
801-
ver))
802-
803782

804783
def _morphed_stc_as_volume(morph, stc, mri_resolution, mri_space, output):
805784
"""Return volume source space as Nifti1Image and/or save to disk."""
806785
assert isinstance(stc, _BaseVolSourceEstimate) # should be guaranteed
807786
if stc._data_ndim == 3:
808787
stc = stc.magnitude()
809-
_check_dep(nibabel='2.1.0', dipy=False)
788+
_check_dep(nibabel='2.1.0')
810789

811790
NiftiImage, NiftiHeader = _triage_output(output)
812791

@@ -926,7 +905,7 @@ def _triage_output(output):
926905

927906
def _interpolate_data(stc, morph, mri_resolution, mri_space, output):
928907
"""Interpolate source estimate data to MRI."""
929-
_check_dep(nibabel='2.1.0', dipy=False)
908+
_check_dep(nibabel='2.1.0')
930909
NiftiImage, NiftiHeader = _triage_output(output)
931910
_validate_type(stc, _BaseVolSourceEstimate, 'stc',
932911
'volume source estimate')
@@ -940,7 +919,7 @@ def _interpolate_data(stc, morph, mri_resolution, mri_space, output):
940919
mri_resolution = (float(mri_resolution),) * 3
941920

942921
if isinstance(mri_resolution, tuple):
943-
_check_dep(nibabel=False, dipy='0.10.1') # nibabel was already checked
922+
_check_dep(dipy='0.10.1') # nibabel was already checked
944923
from dipy.align.reslice import reslice
945924

946925
voxel_size = mri_resolution
@@ -1027,177 +1006,20 @@ def _interpolate_data(stc, morph, mri_resolution, mri_space, output):
10271006
###############################################################################
10281007
# Morph for VolSourceEstimate
10291008

1030-
def _compute_r2(a, b):
1031-
return 100 * (a.ravel() @ b.ravel()) / \
1032-
(np.linalg.norm(a) * np.linalg.norm(b))
1033-
1034-
1035-
def _reslice_normalize(img, zooms):
1036-
from dipy.align.reslice import reslice
1037-
img_zooms = img.header.get_zooms()[:3]
1038-
img_affine = img.affine
1039-
img = _get_img_fdata(img)
1040-
if zooms is not None:
1041-
img, img_affine = reslice(img, img_affine, img_zooms, zooms)
1042-
img /= img.max() # normalize
1043-
return img, img_affine
1044-
1045-
1046-
_ORDERED_STEPS = ('translation', 'rigid', 'affine', 'sdr')
1047-
1048-
1049-
def _validate_niter(niter):
1050-
_validate_type(niter, (dict, list, tuple, None), 'niter')
1051-
niter = _handle_default('transform_niter', niter)
1052-
for key, value in niter.items():
1053-
_check_option('niter key', key, _ORDERED_STEPS)
1054-
_check_option(f'len(niter[{repr(key)}])', len(value), (1, 2, 3))
1055-
return niter
1056-
1057-
1058-
def _compute_morph_sdr(mri_from, mri_to, pipeline, niter, zooms,
1059-
allow_separate_zooms=False):
1009+
def _compute_morph_sdr(moving, static, niter, zooms, prealign=None):
10601010
"""Get a matrix that morphs data from one subject to another."""
1061-
_check_dep(nibabel='2.1.0', dipy='0.10.1')
1062-
from nibabel.spatialimages import SpatialImage
10631011
with np.testing.suppress_warnings():
1064-
from dipy.align import imaffine, imwarp, metrics, transforms
1065-
1066-
#
1067-
# Input validation
1068-
#
1069-
1070-
_validate_type(mri_from, SpatialImage, 'moving image')
1071-
_validate_type(mri_to, SpatialImage, 'static image')
1072-
1073-
# niter, zooms
1074-
types = (list, tuple, 'numeric', None)
1075-
if allow_separate_zooms:
1076-
types = (dict,) + types
1077-
_validate_type(zooms, types, 'zooms')
1078-
zooms = _handle_default('transform_zooms', zooms)
1079-
for key, val in zooms.items():
1080-
_check_option('zooms key', key, _ORDERED_STEPS)
1081-
if val is not None:
1082-
val = tuple(
1083-
float(x) for x in np.array(val, dtype=float).ravel())
1084-
_check_option(f'len(zooms[{repr(key)})', len(val), (1, 3))
1085-
if len(val) == 1:
1086-
val = val * 3
1087-
zooms[key] = val
1088-
niter = _validate_niter(niter)
1089-
1090-
# pipeline
1091-
_validate_type(pipeline, (str, list, tuple), 'pipeline')
1092-
if isinstance(pipeline, str):
1093-
_check_option(
1094-
'pipeline', pipeline, ('all', 'rigids', 'affines'),
1095-
extra='when str')
1096-
if pipeline == 'all':
1097-
pipeline = _ORDERED_STEPS
1098-
elif pipeline == 'rigids':
1099-
pipeline = _ORDERED_STEPS[:_ORDERED_STEPS.index('rigid') + 1]
1100-
else:
1101-
assert pipeline == 'affines'
1102-
pipeline = _ORDERED_STEPS[:_ORDERED_STEPS.index('affine') + 1]
1103-
pipeline = tuple(pipeline)
1104-
for ii, step in enumerate(pipeline):
1105-
name = f'pipeline[{ii}]'
1106-
_validate_type(step, str, name)
1107-
_check_option(name, step, _ORDERED_STEPS)
1108-
ordered_pipeline = tuple(sorted(
1109-
pipeline, key=lambda x: _ORDERED_STEPS.index(x)))
1110-
if pipeline != ordered_pipeline:
1111-
raise ValueError(
1112-
f'Steps in pipeline are out of order, expected {ordered_pipeline} '
1113-
f'but got {pipeline} instead')
1114-
if len(set(pipeline)) != len(pipeline):
1115-
raise ValueError('Steps in pipeline should not be repeated')
1116-
1117-
#
1118-
# Actually do some computations
1119-
#
1120-
if 'sdr' in pipeline:
1121-
kind = 'nonlinear symmetric diffeomorphic registration (SDR)'
1122-
else:
1123-
kind = pipeline[-1].replace('_', ' ') + ' registration'
1124-
logger.info(f'Computing {kind} ...')
1125-
1126-
mri_from_orig = mri_from
1127-
mri_to_orig = mri_to
1128-
last_zooms = -1
1129-
pre_affine = imaffine.AffineMap( # identity transform
1130-
None, mri_to.shape, mri_to.affine, mri_from.shape, mri_to.affine)
1131-
sdr_morph = None # disabled
1132-
mi_metric = imaffine.MutualInformationMetric(nbins=32)
1133-
sigmas = [3.0, 1.0, 0.0]
1134-
factors = [4, 2, 1]
1135-
for step in pipeline:
1136-
assert step in _ORDERED_STEPS
1137-
these_zooms = zooms[step]
1138-
if these_zooms != last_zooms:
1139-
logger.info(f'Reslicing to zooms={these_zooms} for {step} ...')
1140-
mri_from, mri_from_affine = _reslice_normalize(
1141-
mri_from_orig, these_zooms)
1142-
mri_to, mri_to_affine = _reslice_normalize(
1143-
mri_to_orig, these_zooms)
1144-
last_zooms = these_zooms
1145-
# set up Affine Registration
1146-
if step != 'sdr':
1147-
logger.info(f'Optimizing {step}:')
1148-
affreg = imaffine.AffineRegistration(
1149-
mi_metric, level_iters=niter[step],
1150-
sigmas=sigmas, factors=factors)
1151-
if step == 'translation':
1152-
logger.info(' Initializing using center of mass:')
1153-
pre_affine = imaffine.transform_centers_of_mass(
1154-
mri_to, mri_to_affine, mri_from, mri_from_affine)
1155-
dist = np.linalg.norm(pre_affine.affine[:3, 3])
1156-
logger.info(f' Translation: {dist:6.1f} mm')
1157-
trans_inst = transforms.TranslationTransform3D()
1158-
elif step == 'rigid':
1159-
trans_inst = transforms.RigidTransform3D()
1160-
else:
1161-
assert step == 'affine'
1162-
trans_inst = transforms.AffineTransform3D()
1163-
with wrapped_stdout(indent=' ', cull_newlines=True):
1164-
pre_affine = affreg.optimize(
1165-
mri_to, mri_from, trans_inst, None,
1166-
mri_to_affine, mri_from_affine,
1167-
starting_affine=pre_affine.affine)
1168-
del affreg
1169-
else:
1170-
assert step == 'sdr'
1171-
1172-
# Apply the current affine
1173-
mri_from_to = pre_affine.transform(
1174-
mri_from, 'linear', mri_from_affine,
1175-
mri_to.shape, mri_to_affine)
1176-
1177-
# Then run the SDR step (it operates on top of the affine)
1178-
if step == 'sdr':
1179-
logger.info('Optimizing SDR:')
1180-
sdr = imwarp.SymmetricDiffeomorphicRegistration(
1181-
metrics.CCMetric(3), niter['sdr'])
1182-
with wrapped_stdout(indent=' ', cull_newlines=True):
1183-
sdr_morph = sdr.optimize(mri_to, mri_from_to)
1184-
mri_from_to = sdr_morph.transform(mri_from_to)
1185-
1186-
# Report some useful information
1187-
if step in ('translation', 'rigid'):
1188-
dist = np.linalg.norm(pre_affine.affine[:3, 3])
1189-
angle = np.rad2deg(_angle_between_quats(
1190-
np.zeros(3), rot_to_quat(pre_affine.affine[:3, :3])))
1191-
logger.info(f' Translation: {dist:6.1f} mm')
1192-
if step == 'rigid':
1193-
logger.info(f' Rotation: {angle:6.1f}°')
1194-
logger.info(
1195-
f' R²: {_compute_r2(mri_to, mri_from_to):6.1f}%')
1196-
1197-
_debug_img(mri_from, mri_to_affine, 'From-reslice')
1198-
_debug_img(mri_from_to, mri_to_affine, 'From-reslice')
1199-
_debug_img(mri_to, mri_to_affine, 'To-reslice')
1200-
return mri_to_affine, pre_affine, sdr_morph
1012+
from dipy.align import imwarp, metrics
1013+
moving, moving_affine = _reslice_normalize(moving, zooms)
1014+
static, static_affine = _reslice_normalize(static, zooms)
1015+
logger.info(f'Computing nonlinear symmetric diffeomorphic registration '
1016+
'(SDR) ...')
1017+
sdr = imwarp.SymmetricDiffeomorphicRegistration(
1018+
metrics.CCMetric(3), niter)
1019+
with wrapped_stdout(indent=' ', cull_newlines=True):
1020+
sdr_morph = sdr.optimize(static, moving, static_affine, moving_affine,
1021+
prealign=prealign)
1022+
return sdr_morph
12011023

12021024

12031025
def _compute_morph_matrix(subject_from, subject_to, vertices_from, vertices_to,

0 commit comments

Comments
 (0)