[WIP, ENH] Abstracted volume registration #9515

alexrockhill · 2021-06-30T23:23:31Z

This splits up #9503 to just refactor the volume registration for CT/MR without combining it with the existing SDR morph for volume source estimates.

This won't work until mne-tools/mne-misc-data#7 merges.

Needs

Decide on _USE_PREALIGN and equiv in morph.py : one should really work for both (differences are minor)
Add back SourceMorph attribute updates from MRG, ENH: Refactor MRI/CT code #9503

alexrockhill · 2021-07-01T05:07:32Z

After some testing I added

pre_affine, sdr_morph = compute_volume_registration(
    mri_from, mri_to, zooms=zooms,
    niter=dict(translation=niter_affine, rigid=niter_affine,
               affine=niter_affine,
               sdr=niter_sdr if niter_sdr else (1,)))
shape = tuple(sdr_morph.disp_shape)
affine = sdr_morph.disp_grid2world
pre_affine = AffineMap(
    pre_affine, sdr_morph.disp_shape, sdr_morph.disp_grid2world,
    sdr_morph.domain_shape, sdr_morph.domain_grid2world)
if not len(niter_sdr):
    sdr_morph = None

which should give the equivalent outcome to

shape, zooms, affine, pre_affine, sdr_morph = _compute_morph_sdr(
    mri_from, mri_to, niter_affine, niter_sdr, zooms)

but in fact doesn't, here is the stack trace from the tests

(mne) Alexs-MacBook-Pro:mne-python alexrockhill$ pytest mne/tests/test_morph.py
=================================================== test session starts ===================================================
platform darwin -- Python 3.9.4, pytest-6.2.4, py-1.10.0, pluggy-0.13.1
rootdir: /Users/alexrockhill/projects/mne-python, configfile: setup.cfg
plugins: cov-2.12.1, harvest-1.10.3, timeout-1.4.2
collected 596 items                                                                                                       

mne/tests/test_morph.py ........F...F.....FF....................................................................... [ 15%]
................................................................................................................... [ 34%]
................................................................................................................... [ 53%]
................................................................................................................... [ 73%]
................................................................................................................... [ 92%]
.............................................                                                                       [100%]


================================================= slowest 20 test modules =================================================
205.98s total  mne/tests/test_morph.py


======================================================== FAILURES =========================================================
________________________ test_volume_source_morph_round_trip[sample-fsaverage-5.9-6.1-float-False] ________________________
mne/tests/test_morph.py:548: in test_volume_source_morph_round_trip
    assert lower <= mu < upper, f'round-trip distance {mu}'
E   AssertionError: round-trip distance 10.767633310296585
E   assert 10.767633310296585 < 6.1
-------------------------------------------------- Captured stdout call ---------------------------------------------------
    Reading a source space...
    [done]
    1 source spaces read
    Reading a source space...
    [done]
    1 source spaces read
Setting up volume interpolation ...
    21255808/16777216 nonzero values for brain
[done]
Volume source space(s) present...
    Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/mri/brain.mgz as "from" volume
    Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/fsaverage/mri/brain.mgz as "to" volume
Computing registration...
Reslicing to zooms=(7.000000216066837, 7.000000216066837, 7.000000216066837)...
Optimizing translation:
    Optimizing level 0 [max iter: 1]
    Translation:   19.3 mm
    R²:            70.1%
Reslicing to zooms=(7.000000216066837, 7.000000216066837, 7.000000216066837)...
Optimizing rigid:
    Optimizing level 0 [max iter: 1]
    Translation:    0.0 mm
    Rotation:      10.2°
    R²:            90.2%
Reslicing to zooms=(7.000000216066837, 7.000000216066837, 7.000000216066837)...
Optimizing affine:
    Optimizing level 0 [max iter: 1]
    R²:            91.1%
Reslicing to zooms=(7.000000216066837, 7.000000216066837, 7.000000216066837)...
Optimizing sdr:
    R²:            93.6%
[done]
Volume source space(s) present...
    Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/fsaverage/mri/brain.mgz as "from" volume
    Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/mri/brain.mgz as "to" volume
Computing registration...
Reslicing to zooms=(7.000000216066837, 7.000000216066837, 7.000000216066837)...
Optimizing translation:
    Optimizing level 0 [max iter: 1]
    Translation:   19.4 mm
    R²:            70.1%
Reslicing to zooms=(7.000000216066837, 7.000000216066837, 7.000000216066837)...
Optimizing rigid:
    Optimizing level 0 [max iter: 1]
    Translation:    0.0 mm
    Rotation:      18.3°
    R²:            89.8%
Reslicing to zooms=(7.000000216066837, 7.000000216066837, 7.000000216066837)...
Optimizing affine:
    Optimizing level 0 [max iter: 1]
    R²:            94.3%
Reslicing to zooms=(7.000000216066837, 7.000000216066837, 7.000000216066837)...
Optimizing sdr:
    R²:            96.1%
[done]
-------------------------------------------------- Captured stderr call ---------------------------------------------------
100%|██████████| Time : 10/10 [00:00<00:00,   88.46it/s]
100%|██████████| Time : 10/10 [00:00<00:00,   11.11it/s]
_________________________ test_volume_source_morph_round_trip[sample-fsaverage-10-12-float-True] __________________________
mne/tests/test_morph.py:548: in test_volume_source_morph_round_trip
    assert lower <= mu < upper, f'round-trip distance {mu}'
E   AssertionError: round-trip distance 9.288245611270737
E   assert 10 <= 9.288245611270737
-------------------------------------------------- Captured stdout call ---------------------------------------------------
Sphere                : origin at (0.0 0.0 0.0) mm
              radius  : 95.0 mm
grid                  : 5.0 mm
mindist               : 5.0 mm
MRI volume            : /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/mri/aseg.mgz

Reading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/mri/aseg.mgz...

Setting up the sphere...
Surface CM = (   0.0    0.0    0.0) mm
Surface fits inside a sphere with radius   95.0 mm
Surface extent:
    x =  -95.0 ...   95.0 mm
    y =  -95.0 ...   95.0 mm
    z =  -95.0 ...   95.0 mm
Grid extent:
    x = -100.0 ...  100.0 mm
    y = -100.0 ...  100.0 mm
    z = -100.0 ...  100.0 mm
68921 sources before omitting any.
28647 sources after omitting infeasible sources not within 0.0 - 95.0 mm.
24303 sources remaining after excluding the sources outside the surface and less than    5.0 mm inside.
    Selected 12 voxels from 4th-Ventricle
Adjusting the neighborhood info.
Source space : MRI voxel -> MRI (surface RAS)
     0.005000  0.000000  0.000000    -100.00 mm
     0.000000  0.005000  0.000000    -100.00 mm
     0.000000  0.000000  0.005000    -100.00 mm
     0.000000  0.000000  0.000000       1.00
MRI volume : MRI voxel -> MRI (surface RAS)
    -0.003000  0.000000  0.000000     129.00 mm
     0.000000  0.000000  0.003000    -129.00 mm
     0.000000 -0.003000  0.000000     129.00 mm
     0.000000  0.000000  0.000000       1.00
MRI volume : MRI (surface RAS) -> RAS (non-zero origin)
     1.000000 -0.000000 -0.000000      -5.27 mm
    -0.000000  1.000000 -0.000000       9.04 mm
    -0.000000  0.000000  1.000000     -27.29 mm
     0.000000  0.000000  0.000000       1.00
Setting up volume interpolation ...
    1548/636056 nonzero values for 4th-Ventricle
[done]
Sphere                : origin at (0.0 0.0 0.0) mm
              radius  : 95.0 mm
grid                  : 5.0 mm
mindist               : 5.0 mm
MRI volume            : /Users/alexrockhill/mne_data/MNE-testing-data/subjects/fsaverage/mri/aseg.mgz

Reading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/fsaverage/mri/aseg.mgz...

Setting up the sphere...
Surface CM = (   0.0    0.0    0.0) mm
Surface fits inside a sphere with radius   95.0 mm
Surface extent:
    x =  -95.0 ...   95.0 mm
    y =  -95.0 ...   95.0 mm
    z =  -95.0 ...   95.0 mm
Grid extent:
    x = -100.0 ...  100.0 mm
    y = -100.0 ...  100.0 mm
    z = -100.0 ...  100.0 mm
68921 sources before omitting any.
28647 sources after omitting infeasible sources not within 0.0 - 95.0 mm.
24303 sources remaining after excluding the sources outside the surface and less than    5.0 mm inside.
    Selected 16 voxels from 4th-Ventricle
Adjusting the neighborhood info.
Source space : MRI voxel -> MRI (surface RAS)
     0.005000  0.000000  0.000000    -100.00 mm
     0.000000  0.005000  0.000000    -100.00 mm
     0.000000  0.000000  0.005000    -100.00 mm
     0.000000  0.000000  0.000000       1.00
MRI volume : MRI voxel -> MRI (surface RAS)
    -0.001000  0.000000  0.000000     128.00 mm
     0.000000  0.000000  0.001000    -128.00 mm
     0.000000 -0.001000  0.000000     128.00 mm
     0.000000  0.000000  0.000000       1.00
MRI volume : MRI (surface RAS) -> RAS (non-zero origin)
     1.000000  0.000000  0.000000       0.00 mm
     0.000000  1.000000  0.000000       0.00 mm
     0.000000  0.000000  1.000000       0.00 mm
     0.000000  0.000000  0.000000       1.00
Setting up volume interpolation ...
    37304/16777216 nonzero values for 4th-Ventricle
[done]
Volume source space(s) present...
    Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/mri/brain.mgz as "from" volume
    Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/fsaverage/mri/brain.mgz as "to" volume
Computing registration...
Reslicing to zooms=(5.0, 5.0, 5.0)...
Optimizing translation:
    Optimizing level 0 [max iter: 1]
    Translation:   18.6 mm
    R²:            69.9%
Reslicing to zooms=(5.0, 5.0, 5.0)...
Optimizing rigid:
    Optimizing level 0 [max iter: 1]
    Translation:    0.0 mm
    Rotation:      10.5°
    R²:            89.8%
Reslicing to zooms=(5.0, 5.0, 5.0)...
Optimizing affine:
    Optimizing level 0 [max iter: 1]
    R²:            91.6%
Reslicing to zooms=(5.0, 5.0, 5.0)...
Optimizing sdr:
    R²:            92.5%
[done]
Volume source space(s) present...
    Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/fsaverage/mri/brain.mgz as "from" volume
    Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/mri/brain.mgz as "to" volume
Computing registration...
Reslicing to zooms=(5.0, 5.0, 5.0)...
Optimizing translation:
    Optimizing level 0 [max iter: 1]
    Translation:   19.4 mm
    R²:            69.9%
Reslicing to zooms=(5.0, 5.0, 5.0)...
Optimizing rigid:
    Optimizing level 0 [max iter: 1]
    Translation:    0.0 mm
    Rotation:      17.9°
    R²:            89.7%
Reslicing to zooms=(5.0, 5.0, 5.0)...
Optimizing affine:
    Optimizing level 0 [max iter: 1]
    R²:            94.2%
Reslicing to zooms=(5.0, 5.0, 5.0)...
Optimizing sdr:
    R²:            96.0%
[done]
-------------------------------------------------- Captured stderr call ---------------------------------------------------
100%|██████████| Time : 10/10 [00:00<00:00,   64.08it/s]
100%|██████████| Time : 10/10 [00:00<00:00,   24.63it/s]
_______________________________________ test_mixed_source_morph[testing_data-False] _______________________________________
mne/tests/test_morph.py:877: in test_mixed_source_morph
    assert img.astype(bool).sum() == n_want  # correct number of voxels
E   assert 260 == 275
E    +  where 260 = <built-in method sum of numpy.ndarray object at 0x142063030>()
E    +    where <built-in method sum of numpy.ndarray object at 0x142063030> = array([[[[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [Fals...   [[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [False]]]]).sum
E    +      where array([[[[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [Fals...   [[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [False]]]]) = <built-in method astype of numpy.ndarray object at 0x141d04690>(bool)
E    +        where <built-in method astype of numpy.ndarray object at 0x141d04690> = array([[[[0.],\n         [0.],\n         [0.],\n         ...,\n         [0.],\n         [0.],\n         [0.]],\n\n        [[0....     [0.]],\n\n        [[0.],\n         [0.],\n         [0.],\n         ...,\n         [0.],\n         [0.],\n         [0.]]]]).astype
-------------------------------------------------- Captured stdout setup --------------------------------------------------
Setting up the source space with the following parameters:

SUBJECTS_DIR = /Users/alexrockhill/mne_data/MNE-testing-data/subjects
Subject      = sample
Surface      = white
Octahedron subdivision grade 3

>>> 1. Creating the source space...

Doing the octahedral vertex picking...
Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/surf/lh.white...
Mapping lh sample -> oct (3) ...
    Triangle neighbors and vertex normals...
Loading geometry from /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/surf/lh.sphere...
Setting up the triangulation for the decimated surface...
loaded lh.white 66/155407 selected to source space (oct = 3)

Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/surf/rh.white...
Mapping rh sample -> oct (3) ...
    Triangle neighbors and vertex normals...
Loading geometry from /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/surf/rh.sphere...
Setting up the triangulation for the decimated surface...
loaded rh.white 66/156866 selected to source space (oct = 3)

You are now one step closer to computing the gain matrix
Sphere                : origin at (0.0 0.0 0.0) mm
              radius  : 95.0 mm
grid                  : 10.0 mm
mindist               : 5.0 mm
MRI volume            : /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/mri/aseg.mgz

Reading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/mri/aseg.mgz...

Setting up the sphere...
Surface CM = (   0.0    0.0    0.0) mm
Surface fits inside a sphere with radius   95.0 mm
Surface extent:
    x =  -95.0 ...   95.0 mm
    y =  -95.0 ...   95.0 mm
    z =  -95.0 ...   95.0 mm
Grid extent:
    x = -100.0 ...  100.0 mm
    y = -100.0 ...  100.0 mm
    z = -100.0 ...  100.0 mm
9261 sources before omitting any.
3695 sources after omitting infeasible sources not within 0.0 - 95.0 mm.
2969 sources remaining after excluding the sources outside the surface and less than    5.0 mm inside.
    Selected 48 voxels from Left-Cerebellum-Cortex
    Selected 56 voxels from Right-Cerebellum-Cortex
Adjusting the neighborhood info.
Source space : MRI voxel -> MRI (surface RAS)
     0.010000  0.000000  0.000000    -100.00 mm
     0.000000  0.010000  0.000000    -100.00 mm
     0.000000  0.000000  0.010000    -100.00 mm
     0.000000  0.000000  0.000000       1.00
MRI volume : MRI voxel -> MRI (surface RAS)
    -0.003000  0.000000  0.000000     129.00 mm
     0.000000  0.000000  0.003000    -129.00 mm
     0.000000 -0.003000  0.000000     129.00 mm
     0.000000  0.000000  0.000000       1.00
MRI volume : MRI (surface RAS) -> RAS (non-zero origin)
     1.000000 -0.000000 -0.000000      -5.27 mm
    -0.000000  1.000000 -0.000000       9.04 mm
    -0.000000  0.000000  1.000000     -27.29 mm
     0.000000  0.000000  0.000000       1.00
Setting up volume interpolation ...
    36089/636056 nonzero values for Left-Cerebellum-Cortex
    43304/636056 nonzero values for Right-Cerebellum-Cortex
[done]
    Reading a source space...
    [done]
    Reading a source space...
    [done]
    2 source spaces read
Sphere                : origin at (0.0 0.0 0.0) mm
              radius  : 95.0 mm
grid                  : 7.0 mm
mindist               : 5.0 mm
MRI volume            : /Users/alexrockhill/mne_data/MNE-testing-data/subjects/fsaverage/mri/aseg.mgz

Reading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/fsaverage/mri/aseg.mgz...

Setting up the sphere...
Surface CM = (   0.0    0.0    0.0) mm
Surface fits inside a sphere with radius   95.0 mm
Surface extent:
    x =  -95.0 ...   95.0 mm
    y =  -95.0 ...   95.0 mm
    z =  -95.0 ...   95.0 mm
Grid extent:
    x =  -98.0 ...   98.0 mm
    y =  -98.0 ...   98.0 mm
    z =  -98.0 ...   98.0 mm
24389 sources before omitting any.
10443 sources after omitting infeasible sources not within 0.0 - 95.0 mm.
8925 sources remaining after excluding the sources outside the surface and less than    5.0 mm inside.
    Selected 148 voxels from Left-Cerebellum-Cortex
    Selected 127 voxels from Right-Cerebellum-Cortex
Adjusting the neighborhood info.
Source space : MRI voxel -> MRI (surface RAS)
     0.007000  0.000000  0.000000     -98.00 mm
     0.000000  0.007000  0.000000     -98.00 mm
     0.000000  0.000000  0.007000     -98.00 mm
     0.000000  0.000000  0.000000       1.00
MRI volume : MRI voxel -> MRI (surface RAS)
    -0.001000  0.000000  0.000000     128.00 mm
     0.000000  0.000000  0.001000    -128.00 mm
     0.000000 -0.001000  0.000000     128.00 mm
     0.000000  0.000000  0.000000       1.00
MRI volume : MRI (surface RAS) -> RAS (non-zero origin)
     1.000000  0.000000  0.000000       0.00 mm
     0.000000  1.000000  0.000000       0.00 mm
     0.000000  0.000000  1.000000       0.00 mm
     0.000000  0.000000  0.000000       1.00
Volume source space(s) present...
    Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/sample/mri/brain.mgz as "from" volume
    Loading /Users/alexrockhill/mne_data/MNE-testing-data/subjects/fsaverage/mri/brain.mgz as "to" volume
Computing registration...
Reslicing to zooms=(7.0, 7.0, 7.0)...
Optimizing translation:
    Optimizing level 2 [max iter: 1]
    Optimizing level 1 [max iter: 0]
    Optimizing level 0 [max iter: 0]
    Translation:   23.5 mm
    R²:            70.1%
Reslicing to zooms=(7.0, 7.0, 7.0)...
Optimizing rigid:
    Optimizing level 2 [max iter: 1]
    Optimizing level 1 [max iter: 0]
    Optimizing level 0 [max iter: 0]
    Translation:    0.0 mm
    Rotation:       9.8°
    R²:            89.9%
Reslicing to zooms=(7.0, 7.0, 7.0)...
Optimizing affine:
    Optimizing level 2 [max iter: 1]
    Optimizing level 1 [max iter: 0]
    Optimizing level 0 [max iter: 0]
    R²:            91.3%
Reslicing to zooms=(7.0, 7.0, 7.0)...
Optimizing sdr:
    R²:            91.0%
surface source space present ...
Computing morph matrix...
    Left-hemisphere map read.
    Right-hemisphere map read.
    5 smooth iterations done.
    5 smooth iterations done.
[done]
[done]
-------------------------------------------------- Captured stderr call ---------------------------------------------------
100%|██████████| Time : 1/1 [00:00<00:00,   21.37it/s]
_______________________________________ test_mixed_source_morph[testing_data-True] ________________________________________
mne/tests/test_morph.py:877: in test_mixed_source_morph
    assert img.astype(bool).sum() == n_want  # correct number of voxels
E   assert 260 == 275
E    +  where 260 = <built-in method sum of numpy.ndarray object at 0x141f31bd0>()
E    +    where <built-in method sum of numpy.ndarray object at 0x141f31bd0> = array([[[[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [Fals...   [[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [False]]]]).sum
E    +      where array([[[[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [Fals...   [[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [False]]]]) = <built-in method astype of numpy.ndarray object at 0x14206c210>(bool)
E    +        where <built-in method astype of numpy.ndarray object at 0x14206c210> = array([[[[0.],\n         [0.],\n         [0.],\n         ...,\n         [0.],\n         [0.],\n         [0.]],\n\n        [[0....     [0.]],\n\n        [[0.],\n         [0.],\n         [0.],\n         ...,\n         [0.],\n         [0.],\n         [0.]]]]).astype
-------------------------------------------------- Captured stderr call ---------------------------------------------------
100%|██████████| Ori × Time : 3/3 [00:00<00:00,   34.86it/s]
---------------------- generated xml file: /Users/alexrockhill/projects/mne-python/junit-results.xml ----------------------
================================================== slowest 20 durations ===================================================
38.77s call     mne/tests/test_morph.py::test_volume_source_morph_round_trip[sample-sample-0.0-0.1-float-True]
36.00s call     mne/tests/test_morph.py::test_volume_source_morph_round_trip[fsaverage-fsaverage-0.0-0.1-float-False]
30.31s call     mne/tests/test_morph.py::test_volume_source_morph_round_trip[sample-sample-0.0-0.1-complex-False]
28.92s call     mne/tests/test_morph.py::test_volume_source_morph_round_trip[sample-fsaverage-10-12-float-True]
27.24s call     mne/tests/test_morph.py::test_volume_source_morph_basic
25.88s call     mne/tests/test_morph.py::test_volume_source_morph_round_trip[sample-fsaverage-5.9-6.1-float-False]
11.77s setup    mne/tests/test_morph.py::test_mixed_source_morph[testing_data-False]
4.30s call     mne/tests/test_morph.py::test_morph_stc_dense
2.29s call     mne/tests/test_morph.py::test_xhemi_morph
1.39s call     mne/tests/test_morph.py::test_surface_source_morph_round_trip[nearest-0.98-0.99-0-float]
1.32s call     mne/tests/test_morph.py::test_surface_vector_source_morph
1.24s call     mne/tests/test_morph.py::test_surface_source_morph_round_trip[None-0.959-0.963-0-float]
1.05s call     mne/tests/test_morph.py::test_mixed_source_morph[testing_data-True]
0.99s call     mne/tests/test_morph.py::test_mixed_source_morph[testing_data-False]
0.83s call     mne/tests/test_morph.py::test_volume_labels_morph[sl2-88-324-20]
0.75s call     mne/tests/test_morph.py::test_volume_labels_morph[sl0-37-138-8]
0.71s call     mne/tests/test_morph.py::test_volume_labels_morph[sl1-51-204-12]
0.42s call     mne/tests/test_morph.py::test_surface_source_morph_round_trip[3-0.968-0.971-2-complex]
0.36s setup    mne/tests/test_morph.py::test_sourcemorph_consistency
0.09s call     mne/tests/test_morph.py::test_sparse_morph
================================================= short test summary info =================================================
FAILED mne/tests/test_morph.py::test_volume_source_morph_round_trip[sample-fsaverage-5.9-6.1-float-False] - AssertionErr...
FAILED mne/tests/test_morph.py::test_volume_source_morph_round_trip[sample-fsaverage-10-12-float-True] - AssertionError:...
FAILED mne/tests/test_morph.py::test_mixed_source_morph[testing_data-False] - assert 260 == 275
FAILED mne/tests/test_morph.py::test_mixed_source_morph[testing_data-True] - assert 260 == 275
======================================== 4 failed, 592 passed in 222.31s (0:03:42) ========================================

I wrote the matching code by taking examples and matching the output of the shapes and affines to attributes of the sdr_morph.

This suggests that the code that uses affine_registration and applies the SDR with pre_align is functionally different than the code in _compute_morph_sdr which calls each of the component steps individually and then applies the pre_align affine before computing the SDR.

I think the accuracy of the affine registration and SDR pipeline in this PR is superior to _compute_morph_sdr because I've tested it thoroughly and it works well for aligning a CT to an MR and for morphing a subject's brain to the fsaverage brain. The _compute_morph_sdr has several mistakes that point to a lack of attention to detail such as passing zooms and then returning it unaltered and the general unnecessary complexity and lack of succinctness that make me hesitant to look into it and see if it's better, which seems unlikely. I think it would be better to switch to the code in this PR in a subsequent PR that fixes the four failing tests and cleans things up a bit. The code in this PR is cleaner, more succinct and works so I think it would be better to go with it.

larsoner · 2021-07-01T12:32:07Z

I think the accuracy of the affine registration and SDR pipeline in this PR is superior to _compute_morph_sdr because I've tested it thoroughly and it works well for aligning a CT to an MR and for morphing a subject's brain to the fsaverage brain.

This makes sense, and I agree this is strong evidence that it works!

The _compute_morph_sdr has several mistakes that point to a lack of attention to detail such as passing zooms and then returning it unaltered and the general unnecessary complexity and lack of succinctness that make me hesitant to look into it and see if it's better, which seems unlikely.

Here I think you sell the current code a bit short. That code has been around for a long time, and had probably half a dozen bugs fixed. This has led some some cruft, as you notice, but it also has led to new unit tests each time to check to make sure things are working properly. Hence the failing tests -- especially the round-trip distance one -- suggest that something isn't as good in the current implementation.

To me there are two issues / potential differences between this (and the custom code in main) and #9503 -- the quality/methods of the affine (MR-CT) alignment, and then how the affine is applied for the SDR morph calculation (i.e., whether it's given the prealign, or affine-transformed-then-computed). I think it's important we somewhat separate these issues when thinking about this.

From what you're saying, it sounds to me like you have two objections to #9503:

The MR-CT alignment is off (this is rigid, does not involve SDR).
The SDR morph (which builds on top of the affine pipeline) works okay in that example, but are not confident in it compared to the prealign case, which you have tested extensively.

Is that right? I agree we should improve the rigid/affine alignment (1), and if this PR does that then it seems to be on the right track. But I'm not sure whether or not prealign will be preferable to the apply-then-compute SDR morph (2). I think it remains to be shown that it works better for the SDR step, especially since neither this nor #9503 allow using the improved affine code here and non-prealigned SDR. I push this point because I suspect it might be causing at least some of the test failures--in particular I seem to remember the round-trip voxel count 276/260 being due to prealign rather than apply-affine-and-compute-SDR.

My suggestion to get to the bottom of it is to look at why the tests fail, and see if they're fixed by not using apply-affine-then-compute-SDR rather than prealign-via-compute -- and then decide whether or not they actually show something is working functionally better or not (they really should since they were introduced to fix bugs, but who knows!).

* upstream/main: MAINT: Speed up CIs (mne-tools#9518) Update tools (mne-tools#9517) FIX: Fix simulate_evoked and apply_forward (mne-tools#9513) FIX: Mayavi (mne-tools#9512)

alexrockhill · 2021-07-01T16:45:52Z

Sounds good, thanks for the feedback, I'll take a look into the failing tests.

larsoner · 2021-07-01T16:46:19Z

Actually I think your snippet above just needs a small tweak, I think I almost have it working. Give me 10 minutes...

alexrockhill · 2021-07-01T16:48:03Z

Actually I think your snippet above just needs a small tweak, I think I almost have it working. Give me 10 minutes...

Hey yay teamwork, take your time!

larsoner · 2021-07-01T19:29:14Z

Whenever we're settled here, I want to port the remaining morph changes over that give finer-grained control to niter. So I'm going to set this to WIP

alexrockhill · 2021-07-01T21:27:37Z

Decide on _USE_PREALIGN and equiv in morph.py : one should really work for both (differences are minor)

So I think it's clear that although the prealign has a bit weirder behavior, it is the only practical option for warping an indexed elec_image. This version that I just pushed has a failing tests that depends on the SDR with aligned input:

--------------------------------------------------------------- Captured stderr call ----------------------------------------------------------------
100%|██████████| Time : 1/1 [00:00<00:00,   21.51it/s]
____________________________________________________ test_mixed_source_morph[testing_data-True] _____________________________________________________
mne/tests/test_morph.py:877: in test_mixed_source_morph
    assert img.astype(bool).sum() == n_want  # correct number of voxels
E   assert 266 == 275
E    +  where 266 = <built-in method sum of numpy.ndarray object at 0x13ced5390>()
E    +    where <built-in method sum of numpy.ndarray object at 0x13ced5390> = array([[[[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [Fals...   [[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [False]]]]).sum
E    +      where array([[[[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [Fals...   [[False],\n         [False],\n         [False],\n         ...,\n         [False],\n         [False],\n         [False]]]]) = <built-in method astype of numpy.ndarray object at 0x13ced50f0>(bool)
E    +        where <built-in method astype of numpy.ndarray object at 0x13ced50f0> = array([[[[0.],\n         [0.],\n         [0.],\n         ...,\n         [0.],\n         [0.],\n         [0.]],\n\n        [[0....     [0.]],\n\n        [[0.],\n         [0.],\n         [0.],\n         ...,\n         [0.],\n         [0.],\n         [0.]]]]).astype
--------------------------------------------------------------- Captured stderr call ----------------------------------------------------------------
100%|██████████| Ori × Time : 3/3 [00:00<00:00,   32.58it/s]

that I'm not super inclined to look into fixing, but since the differences are minor maybe you could do that last fix @larsoner. Alternately, I'm about to push a version that splits the SDR morph in morph.py back and that can be a temporary fix until another PR happens to unify the two morphs. Feel free to revert the last commit if you want to try and fit them together.

…ll pass

alexrockhill · 2021-07-02T16:47:54Z

Just running the tutorial on the new new PR but looks like it works

alexrockhill added 2 commits June 30, 2021 16:21

abstracted volume registration

d7ebf55

fix doc

59a239a

Merge remote-tracking branch 'upstream/main' into refactor

d79c929

* upstream/main: MAINT: Speed up CIs (mne-tools#9518) Update tools (mne-tools#9517) FIX: Fix simulate_evoked and apply_forward (mne-tools#9513) FIX: Mayavi (mne-tools#9512)

ENH: Unify approaches

adb72ad

larsoner force-pushed the refactor branch from bea6109 to adb72ad Compare July 1, 2021 19:25

larsoner changed the title ~~[MRG, ENH] Abstracted volume registration~~ [WIP, ENH] Abstracted volume registration Jul 1, 2021

larsoner mentioned this pull request Jul 1, 2021

MRG, ENH: Refactor MRI/CT code #9503

Closed

7 tasks

larsoner and others added 5 commits July 1, 2021 16:07

FIX: Wrong order

07c6b0a

FIX: Justified

f771fb8

get tutorial working

653dd3b

prealign version, working all but two tests

fc9b9c3

typo

bda2f05

alexrockhill added 5 commits July 1, 2021 15:02

merge this version with minor changes, separate PR to combine

27a15e9

last bit of cleanup

9a8087d

arggg struggling with tests, had to loosen/change two distances but a…

672ac70

…ll pass

one more very small test loosen

9d67ef4

one more test loosening

729ae8a

larsoner mentioned this pull request Jul 2, 2021

MRG, ENH: Abstracted volume registration again #9521

Merged

alexrockhill closed this Jul 2, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

[WIP, ENH] Abstracted volume registration #9515

[WIP, ENH] Abstracted volume registration #9515

Uh oh!

alexrockhill commented Jun 30, 2021 •

edited by larsoner

Loading

Uh oh!

alexrockhill commented Jul 1, 2021 •

edited

Loading

Uh oh!

larsoner commented Jul 1, 2021

Uh oh!

alexrockhill commented Jul 1, 2021

Uh oh!

larsoner commented Jul 1, 2021

Uh oh!

alexrockhill commented Jul 1, 2021

Uh oh!

larsoner commented Jul 1, 2021

Uh oh!

alexrockhill commented Jul 1, 2021 •

edited

Loading

Uh oh!

alexrockhill commented Jul 2, 2021

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

[WIP, ENH] Abstracted volume registration #9515

[WIP, ENH] Abstracted volume registration #9515

Uh oh!

Conversation

alexrockhill commented Jun 30, 2021 • edited by larsoner Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

alexrockhill commented Jul 1, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

larsoner commented Jul 1, 2021

Uh oh!

alexrockhill commented Jul 1, 2021

Uh oh!

larsoner commented Jul 1, 2021

Uh oh!

alexrockhill commented Jul 1, 2021

Uh oh!

larsoner commented Jul 1, 2021

Uh oh!

alexrockhill commented Jul 1, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

alexrockhill commented Jul 2, 2021

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

alexrockhill commented Jun 30, 2021 •

edited by larsoner

Loading

alexrockhill commented Jul 1, 2021 •

edited

Loading

alexrockhill commented Jul 1, 2021 •

edited

Loading