MRG, ENH: Refactor MRI/CT code #9503

larsoner · 2021-06-24T18:18:52Z

@alexrockhill can you take over the first two steps below:

Fix example so that all electrodes get assigned locations (currently the np.where can be length zero, which is bad)
Make sure example is still structurally okay

Then I can do for the functions I added:

Reduce footprint from 71.64 sec 2537.8 MB
Add tests
Speed up tests
Add imposed translation, ensure translation is found in each pipeline
Update SourceMorph object attributes niter_* -> niter dict, add pipeline

larsoner · 2021-06-24T18:25:43Z

Oh wow 71.64 sec 2537.8 MB that's too much memory. I'll have to figure out how to reduce that...

larsoner · 2021-06-24T18:35:18Z

Just pushed a commit to make CircleCI complete without error (my hack workaround for having none at the right value was wrong), feel free to push right over it @alexrockhill if you have a correct fix.

alexrockhill · 2021-06-24T18:36:51Z

Hmmm, yeah I just tried and it doesn't work either. I think that's because the rigid registration can't be zoomed... Do you want to split it up into two steps? mne.morph.compute_volume_translation and mne.morph.compute_volume_alignment... I don't think warp is that accurate either, you're not changing the scaling, just rotating and translating

alexrockhill · 2021-06-24T18:39:30Z

The reason np.where is empty is because the elec_image that got morphed was off because the alignment was off.

larsoner · 2021-06-24T18:41:52Z

I don't think warp is that accurate either, you're not changing the scaling, just rotating and translating

You must mean for the CT-MR warp step in the subject's own space. I'll I add a rigid=False param that skips the affine step that only does trans+rot.

I don't think warp is that accurate either, you're not changing the scaling, just rotating and translating

Is this still in reference to the CT-MR warp, or the MR-template warp?

alexrockhill · 2021-06-24T18:43:29Z

I don't think warp is that accurate either, you're not changing the scaling, just rotating and translating

You must mean for the CT-MR warp step in the subject's own space. I'll I add a rigid=False param that skips the affine step that only does trans+rot.

I don't think warp is that accurate either, you're not changing the scaling, just rotating and translating

Is this still in reference to the CT-MR warp, or the MR-template warp?

In reference to the MR-CT registration which I don't think qualifies as a warp and so I think is misleading.

alexrockhill · 2021-06-24T18:44:34Z

I think it should be split into two separate registration and warp steps, I can take a pass if you'd like

larsoner · 2021-06-24T18:50:45Z

I think it should be split into two separate registration and warp steps, I can take a pass if you'd like

To me it's all registration/warping/coreg, just with different levels of granularity (translation is one, rotation is another, affine-but-still-linear / shear / scaling is another, SDR is another). This is why I'm inclined toward a single function to do it...

larsoner · 2021-06-24T18:51:39Z

I made the first step rigid and it still has invalid electrodes

alexrockhill · 2021-06-24T18:54:29Z

I made the first step rigid and it still has invalid electrodes

If it has the same zoom=5 then it won't align properly, the CT to MR rigid can't be zoomed.

alexrockhill · 2021-06-24T18:55:49Z

I think it should be split into two separate registration and warp steps, I can take a pass if you'd like

To me it's all registration/warping/coreg, just with different levels of granularity (translation is one, rotation is another, affine-but-still-linear / shear / scaling is another, SDR is another). This is why I'm inclined toward a single function to do it...

I just think when you're passing a CT into a function to align with an MR, you don't want it to be "warped", that would be bad, just registered from my best guess at a naive user's perspective

alexrockhill · 2021-06-24T18:59:10Z

Maybe a good solution is to call it mne.transforms.register with a sdr_warp=True keyword argument to do the SDR as well all in one call. I think even having under in the morph module is misleading since you don't want to "morph" the CT either.

That's just my perspective though, maybe I'm being unreasonable and others are perfectly fine having rigid registrations under morph and warp.

larsoner · 2021-06-24T19:00:31Z

Maybe under mne.source_space is best actually, that's where a bunch of MRI wrangling things live

larsoner · 2021-06-24T19:02:06Z

... I think you're right that mne.transforms is probably best, I'll push

larsoner · 2021-06-24T19:27:20Z

Okay all electrodes get assigned now, feel free to look at the outputs to see if they make sense, and also the naming/structure etc.

alexrockhill · 2021-06-24T19:39:06Z

Looks great!

I don't have write access but I would suggest having zooms be consistent with niter and either have zooms_trans, zooms_affine and zooms_sdr or have one niter dictionary. I would break up zooms_trans and zooms_affine as that will speed things up for the CT registration; it's nice to zoom the translation and not the rigid affine for how I would actually do it for a real analysis.

larsoner · 2021-06-24T19:42:30Z

. I would break up zooms_trans and zooms_affine as that will speed things up for the CT registration;

This suggests we might actually want the entire set zooms_translation, zooms_rigid, zooms_affine, zooms_sdr. I'm inclined to make both niter and zooms dicts for this reason, we'll end up with a ton of parameters otherwise. And maybe the logic should be that if a given zooms level is not provided, the zooms of the previous level are used. If zooms['whatever'] = None it means "use native". So for example zooms=dict(translation=5, rigid=None) would mean "use 5 for translation, then native from rigid onward"

alexrockhill · 2021-06-24T19:42:43Z

I think also a keyword warp for whether to compute the SDR would be really helpful, to me niters_sdr=() is confusing

larsoner · 2021-06-24T19:44:11Z

I think also a keyword warp for whether to compute the SDR would be really helpful, to me niters_sdr=() is confusing

How about niter=dict(affine=None, sdr=()) would mean "only do rigid" (i.e., use None as an alias for () meaning, "do not do it")?

alexrockhill · 2021-06-24T19:44:53Z

. I would break up zooms_trans and zooms_affine as that will speed things up for the CT registration;

This suggests we might actually want the entire set zooms_translation, zooms_rigid, zooms_affine, zooms_sdr. I'm inclined to make both niter and zooms dicts for this reason, we'll end up with a ton of parameters otherwise. And maybe the logic should be that if a given zooms level is not provided, the zooms of the previous level are used. If zooms['whatever'] = None it means "use native". So for example zooms=dict(translation=5, rigid=None) would mean "use 5 for translation, then native from rigid onward"

That sounds great except that I wouldn't want it to carry on at a lower resolution. If the previous level used some downsampling, that doesn't mean the next level has to, in fact in the case of the CT to MR, you want to zoom and then not. On None I would just use the full resolution.

alexrockhill · 2021-06-24T19:47:06Z

I think also a keyword warp for whether to compute the SDR would be really helpful, to me niters_sdr=() is confusing

How about niter=dict(affine=None, sdr=()) would mean "only do rigid" (i.e., use None as an alias for () meaning, "do not do it")?

I think that's getting really confusing to me. Maybe a steps list would be good, with ['trans', 'rigid'] for CT to MR and ['trans', 'rigid', 'affine', 'sdr'] for the T1 to fsaverage.

In Dipy it's pipeline so I guess that would be appropriate too.

Should be a tuple if the default is ('trans', 'rigid', 'affine', 'sdr') so that it's immutable also.

larsoner · 2021-06-24T19:52:04Z

I think that's getting really confusing to me. Maybe a steps list would be good, with ['trans', 'rigid'] for CT to MR and ['trans', 'rigid', 'affine', 'sdr'] for the T1 to fsaverage.

That seems like a nice way to do it!

That sounds great except that I wouldn't want it to carry on at a lower resolution. If the previous level used some downsampling, that doesn't mean the next level has to, in fact in the case of the CT to MR, you want to zoom and then not. On None I would just use the full resolution.

So the logic would be "native unless specified"? I can live with this, especially if we alias zooms=5 to mean "use 5 for all steps". Even in the example here we use 6 then 4 for subject-to-fsaverage, I'll probably switch it to zooms=5 in that case (and the results are still fine it seems).

alexrockhill · 2021-06-24T19:56:34Z

I think that's getting really confusing to me. Maybe a steps list would be good, with ['trans', 'rigid'] for CT to MR and ['trans', 'rigid', 'affine', 'sdr'] for the T1 to fsaverage.

That seems like a nice way to do it!

That sounds great except that I wouldn't want it to carry on at a lower resolution. If the previous level used some downsampling, that doesn't mean the next level has to, in fact in the case of the CT to MR, you want to zoom and then not. On None I would just use the full resolution.

So the logic would be "native unless specified"? I can live with this, especially if we alias zooms=5 to mean "use 5 for all steps". Even in the example here we use 6 then 4 for subject-to-fsaverage, I'll probably switch it to zooms=5 in that case (and the results are still fine it seems).

Sounds perfect, totally on the same page

larsoner · 2021-06-24T19:57:10Z

Okay I'll get to it tomorrow, or if you want to give it a shot you can

alexrockhill · 2021-06-25T06:12:42Z

I didn't finish and I need to fix the commit history but I pushed these in case you want to finish it tomorrow morning before I'm up, but I can also fix and finish tomorrow just a bit later, I'm on PST.

doc/changes/latest.inc

agramfort · 2021-06-25T07:16:59Z

mne/transforms.py

+
+
+@verbose
+def compute_volume_registration(moving, static, steps='warp', zooms=None,


I would call moving_img and static_img to suggest it's a nibabel Image

I had it that way at first but thought moving and static alone were clearer because it matches what dipy calls it in affine_registration:

https://dipy.org/documentation/1.4.1./reference/dipy.align/#affine-registration

So here I think being consistent with their naming is worth it. (Plus the name is shorter.) But I don't feel strongly so if you aren't convinced by the consistency argument I'll change it

mne/viz/topomap.py

mne/transforms.py

This reverts commit 6cd2500.

This reverts commit 22f3f14.

This reverts commit 07c07e3.

This reverts commit 288cd91.

This reverts commit 97c43cd.

This reverts commit be9154a.

This reverts commit cb132a6.

fix core fix affine map working version for tutorial, not tests fix test

This reverts commit 50e8107.

larsoner · 2021-06-30T17:20:22Z

Okay @alexrockhill feel free to have a look and see if you're happy with the code, here is the updated example:

https://30168-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/clinical/10_ieeg_localize.html

larsoner · 2021-06-30T18:24:59Z

Windows is having a bad day because of the VTK 9.0.2 release, can be ignored here as I'm working on it in #9512

alexrockhill · 2021-06-30T19:07:12Z

Okay @alexrockhill feel free to have a look and see if you're happy with the code, here is the updated example:

https://30168-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/clinical/10_ieeg_localize.html

Still looks off... I would rather have a maximally accurate example that has a lot of code that's not abstracted than one that doesn't work quite right. I think the existing morph functions may not be compatible with the what is on main now for the tutorial. I think first priority if you want to go this route of refactoring is getting a version that works and then maybe in a followup we can get the SDR morph from the volume source estimate working using the same code. I really would rather not spend a ton of time debugging the volume source estimate code so I think it's best to part ways. I'll take a pass now.

larsoner · 2021-06-30T19:11:11Z

Still looks off...

Can you clarify how/where? I'm assuming you're talking about subject->fsaverage and not CT-MRI. For "to fsaverage", those temporal electrodes are now nicely on/around the actual lobe, which is what I was always worried about, so to me it looks okay. Is there an image above somewhere (or in some artifact) that looks better?

I'll take a pass now.

Can you see if just zooms=None in the SDR pass makes things look okay to you? If so I can go back to saving the warped MRIs.

I really would rather not spend a ton of time debugging the volume source estimate code so I think it's best to part ways.

So maybe as an intermediate, if zooms=None does not work, we can keep the MR-CT alignment that is here currently since it works, then put back your custom dipy code for the subject->fsaverage alignment and save those warped electrode and brain images?

alexrockhill · 2021-06-30T19:12:07Z

Still looks off...

Can you clarify how/where? I'm assuming you're talking about subject->fsaverage and not CT-MRI. For "to fsaverage", those temporal electrodes are now nicely on/around the actual lobe, which is what I was always worried about, so to me it looks okay. Is there an image above somewhere (or in some artifact) that looks better?

I'll take a pass now.

Can you see if just zooms=None in the SDR pass makes things look okay to you? If so I can go back to saving the warped MRIs.

I really would rather not spend a ton of time debugging the volume source estimate code so I think it's best to part ways.

So maybe as an intermediate, if zooms=None does not work, we can keep the MR-CT alignment that is here currently since it works, then put back your custom dipy code for the subject->fsaverage alignment and save those warped electrode and brain images?

The CT-MR is way off still, that's the main thing. The SDR looks pretty good.

larsoner · 2021-06-30T19:18:27Z

The CT-MR is way off still, that's the main thing.

So this is what we have now:

And this is what we have in main:

I think it only differs by 1.2mm and 1.9 degrees from what you had before. I have a hard time seeing meaningful differences here, and if I didn't know which is which I'm not sure which I'd say was better, but I am no expert here!

As much as I'm inclined to keep it, indeed we shouldn't keep something if you think it is "way off" even if I can't see it. I thought you had said before that the MR-CT was now okay, but since it's not I agree perhaps the easiest thing to do is just to revert to the old affine and commented out code for the MR-CT step, and revisit sometime later how to use our standard pipeline to get an affine of equivalent quality.

alexrockhill · 2021-06-30T19:22:41Z

The CT-MR is way off still, that's the main thing.

So this is what we have now:

And this is what we have in main:

I think it only differs by 1.2mm and 1.9 degrees from what you had before. I have a hard time seeing meaningful differences here, and if I didn't know which is which I'm not sure which I'd say was better, but I am no expert here!

As much as I'm inclined to keep it, indeed we shouldn't keep something if you think it is "way off" even if I can't see it. I thought you had said before that the MR-CT was now okay, but since it's not I agree perhaps the easiest thing to do is just to revert to the old affine and commented out code for the MR-CT step, and revisit sometime later how to use our standard pipeline to get an affine of equivalent quality.

Way off is a relative term and I think you have to spend some time aligning these by hand to appreciate when the fit is right and when it's not. 1.2 mm and 1.9 degrees isn't a whole lot but those are clearly how much it's wrong and to me, that's too much. You can align by eye closer than that, the whole point of this is that it should be as precise as possible with the alignment, better than what can be done by human eye. I could fix the top image to be much closer by hand.

alexrockhill · 2021-06-30T19:24:49Z

Keep in mind that you're only seeing a bit of the alignment with the three slices and in a viewer like freeview it would be much more obvious how much the alignment is actually off and I would consider the first alignment to be poor and need fixing if someone had done it by hand and I was reviewing it and the second alignment would be good.

larsoner · 2021-06-30T19:27:45Z

tutorials/clinical/10_ieeg_localize.py

-alignment_affine = np.array([
-    [0.99235816, -0.03412124, 0.11857915, -133.22262329],
-    [0.04601133, 0.99402046, -0.09902669, -97.64542095],
-    [-0.11449119, 0.10372593, 0.98799428, -84.39915646],


@alexrockhill okay if I just revert the text removal here, then we're good to merge?

That's the correct alignment but the commented out code should return that matrix and it doesn't right?

By revert the text I mean both this affine and the old commented code, which I assume is what you used to generate it

alexrockhill · 2021-06-30T23:26:40Z

Ok I just pushed a different pull request to split this up. First, I think a version that works with the tutorial should be merged which is the PR I just pushed. Then we can try and refactor the volume source estimate code to use the new transforms. I think combining the two steps is causing too much of an issue.

larsoner · 2021-07-01T19:30:45Z

Closing for #9515

agramfort reviewed Jun 25, 2021

View reviewed changes

larsoner commented Jun 25, 2021

View reviewed changes

mne/transforms.py Outdated Show resolved Hide resolved

larsoner added 3 commits June 25, 2021 09:20

ENH: Refactor MRI/CT code

2fa62e2

FIX: Allow to complete

5b5173c

FIX: Allow rigid

124fc3f

larsoner and others added 12 commits June 29, 2021 15:51

Revert "cut out multiple linear interpolations"

93cf9a6

This reverts commit 6cd2500.

Revert "put back linear interpolation on all steps"

6baa5c9

This reverts commit 22f3f14.

Revert "affine registration fix"

48a13fd

This reverts commit 07c07e3.

Revert "wrong affine in apply"

45bc345

This reverts commit 288cd91.

Revert "forgot blank line"

ba33563

This reverts commit 97c43cd.

Revert "flake"

d472500

This reverts commit be9154a.

Revert "working version, lot's of reverting"

98dc4d3

This reverts commit cb132a6.

FIX: Example

0a6413f

working version for tutorial, not tests

50e8107

fix core fix affine map working version for tutorial, not tests fix test

WIP: Uncomment [ci skip]

2be42d3

Revert "working version for tutorial, not tests"

d9fc2e2

This reverts commit 50e8107.

FIX: Fix affines

aebfbc9

larsoner changed the title ~~WIP, ENH: Refactor MRI/CT code~~ MRG, ENH: Refactor MRI/CT code Jun 30, 2021

larsoner mentioned this pull request Jun 30, 2021

Add sphinx-gallery Gallery to the website vispy/vispy#2122

Closed

larsoner commented Jun 30, 2021

View reviewed changes

alexrockhill mentioned this pull request Jun 30, 2021

[WIP, ENH] Abstracted volume registration #9515

Closed

2 tasks

larsoner mentioned this pull request Jul 1, 2021

MRG, MAINT: Speed up CIs #9518

Merged

larsoner closed this Jul 1, 2021

larsoner mentioned this pull request Jul 2, 2021

MRG, ENH: Abstracted volume registration again #9521

Merged

larsoner deleted the tut branch February 1, 2022 16:12



		@verbose
		def compute_volume_registration(moving, static, steps='warp', zooms=None,

Uh oh!

MRG, ENH: Refactor MRI/CT code #9503

MRG, ENH: Refactor MRI/CT code #9503

Conversation

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

Uh oh!

Uh oh!

larsoner commented Jun 24, 2021

Uh oh!

larsoner commented Jun 24, 2021

Uh oh!

alexrockhill commented Jun 24, 2021

Uh oh!

alexrockhill commented Jun 24, 2021

Uh oh!

larsoner commented Jun 24, 2021

Uh oh!

alexrockhill commented Jun 24, 2021

Uh oh!

alexrockhill commented Jun 24, 2021

Uh oh!

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

Uh oh!

Uh oh!

larsoner commented Jun 24, 2021

Uh oh!

alexrockhill commented Jun 24, 2021

Uh oh!

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

Uh oh!

Uh oh!

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

Uh oh!

Uh oh!

larsoner commented Jun 24, 2021

Uh oh!

larsoner commented Jun 24, 2021

Uh oh!

larsoner commented Jun 24, 2021

Uh oh!

alexrockhill commented Jun 24, 2021

Uh oh!

larsoner commented Jun 24, 2021

Uh oh!

alexrockhill commented Jun 24, 2021

Uh oh!

larsoner commented Jun 24, 2021

Uh oh!

alexrockhill commented Jun 24, 2021

Uh oh!

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

Uh oh!

Uh oh!

larsoner commented Jun 24, 2021

Uh oh!

alexrockhill commented Jun 24, 2021

Uh oh!

larsoner commented Jun 24, 2021

Uh oh!

alexrockhill commented Jun 25, 2021

Uh oh!

Uh oh!

agramfort Jun 25, 2021

Choose a reason for hiding this comment

Uh oh!

larsoner Jun 25, 2021

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

larsoner commented Jun 30, 2021

Uh oh!

larsoner commented Jun 30, 2021

Uh oh!

alexrockhill commented Jun 30, 2021

Uh oh!

larsoner commented Jun 30, 2021

Uh oh!

alexrockhill commented Jun 30, 2021

Uh oh!

larsoner commented Jun 24, 2021 •

edited

Loading

larsoner commented Jun 24, 2021 •

edited

Loading

alexrockhill commented Jun 24, 2021 •

edited

Loading

alexrockhill commented Jun 24, 2021 •

edited

Loading

alexrockhill commented Jun 24, 2021 •

edited

Loading