Improve moco for interleaved data

[valosekj]

Added -interleaved option to sct_dmri_moco.py, if it is selected, newly created function in moco.py for splitting input data to even and odd datasets is called. This function splits data and calls moco_wrapper function two times, then merge datasets back together.
This current solution is created to influence original code as least as possible and needs improvements.

I tried current workaround on testing data sub-issue2637 with commands from http://forum.spinalcordmri.org/t/sct-dmri-moco-doesnt-seem-to-correct-motion/358/6 (average DWI, segment SC, create mask, run moco).
Compare to results from original moco, I can see improvements in sagittal plane in slice-by-slice misalignment.

Fixes #2637

[kousu]

Heads up, this might run into #2642. It doesn't edit sct_dmri_moco directly but it does change the underlying image registration algorithm; I had to update the snapshot tests to get CI to pass -- not by much, but by enough. When that gets merged you should rebase on top of it ((including reinstalling sct_download_data -d binaries_linux -o bin/)) early to catch any issues.

[valosekj]

@valosekj wasn't the point of #2684 to be able to develop a specific test to quantify the improvements of your moco-interleaved feature?

You are right! I meant by "Tested on sub-glen." that all parts of code work and script finish without errors. I'll test/validate performance of this PR on some systemic data asap.

[jcohenadad]

@valosekj wasn't the point of #2684 to be able to develop a specific test to quantify the improvements of your moco-interleaved feature?

You are right! I meant by "Tested on sub-glen." that all parts of code work and script finish without errors. I'll test/validate performance of this PR on some systemic data asap.

ok let me know when the PR is at its final stage for me to review

[valosekj]

I have tested this PR on a synthetic 4D dummy interleaved data generated by create_test_data.py.
I run both moco (sct_dmri_moco -interleaved 0) and moco_interleaved (sct_dmri_moco -interleaved 1) on this data. I can see that generated moco_params.tsv files differ between original moco and moco_interleaved. But visually in fsleyes, I can see only very slight (or mostly no) difference between raw image and images after moco or moco_interleaved.
Therefore, I am now a little bit unsure how exactly validate efficiency of the interleaved approach.

Code how to generate synthetic interleaved 4D data:

from spinalcordtoolbox.testing.create_test_data import dummy_segmentation_4d
DEBUG = True  # Set to True to save images
seg_4D = dummy_segmentation_4d(vol_num=10, create_bvecs=True, size_arr=(64, 64, 49), pixdim=(1,1,1),radius_RL=13.0, radius_AP=5.0, shape='ellipse', interleaved=True, debug=DEBUG)

dummy_4D_interleaved.zip

[jcohenadad]

But visually in fsleyes, I can see only very slight (or mostly no) difference between raw image and images after moco or moco_interleaved.
Therefore, I am now a little bit unsure how exactly validate efficiency of the interleaved approach.

what if you try to increase the amount of motion in your simulation? how much motion is there now?

[jcohenadad]

i've looked at the simulated data uploaded in #2664 (comment), motion is very small (sub-mm) in most of the middle slices:

so i suggest you increase motion to be able to better validate the performance of your algorithm @valosekj

also: i suggest you upload gif anim-- more convenient that downloading zip, opening fsleyes, etc.

[valosekj]

But visually in fsleyes, I can see only very slight (or mostly no) difference between raw image and images after moco or moco_interleaved.
Therefore, I am now a little bit unsure how exactly validate efficiency of the interleaved approach.

what if you try to increase the amount of motion in your simulation? how much motion is there now?

Generated image has 1mm isotropic resolution (or should have, because, to be precise - pixdim flag is set to (1,1,1):

https://github.com/neuropoly/spinalcordtoolbox/blob/67c679deadd8af0c2eb7445041b0c61628d8bbfa/spinalcordtoolbox/testing/create_test_data.py#L217
but fsleyes states at pixdim "[unknown units]"):

Slice-wise shift is in range between 0.5-3mm (or again - unknown units?). Unfortunately now set as fixed value (I should change it to be input arg):
https://github.com/neuropoly/spinalcordtoolbox/blob/67c679deadd8af0c2eb7445041b0c61628d8bbfa/spinalcordtoolbox/testing/create_test_data.py#L232

(for context; middle slices have lower shift that outer ones due to polynomial function - see this image)

I'll increase motion and will update gif, thanks 👍

[jcohenadad]

middle slices have lower shift that outer ones due to polynomial function

that's not necessarily the case, e.g. with 2nd-order you could have largest shift in the middle (depending on how you define you constant 0-term offset)

[valosekj]

middle slices have lower shift that outer ones due to polynomial function

that's not necessarily the case, e.g. with 2nd-order you could have largest shift in the middle (depending on how you define you constant 0-term offset)

hm, currently, the order of polynomial is equal to number of slices (nz):
https://github.com/neuropoly/spinalcordtoolbox/blob/67c679deadd8af0c2eb7445041b0c61628d8bbfa/spinalcordtoolbox/testing/create_test_data.py#L160

Lower degree results into worse fit (i.e. slicewise regularized fit (orange) does not properly fit/simulate interleaved acquisition mode (blue dots)):

[jcohenadad]

@valosekj what exactly are the graphs #2664 (comment) supposed to represent?

if the blue dots are the pixel shift in the Y direction across slices, for both odd and even slices, then why are you trying to fit both odd and even slices at the same time? shouldn't the simulated data be constructed by fitting one polynomial function to the odd slices, and another one to the even slices?

[valosekj]

if the blue dots are the pixel shift in the Y direction across slices, for both odd and even slices

Yes, blue dots represent the pixel shift in A-P (Y) direction across slices (S-I direction (Z)) for one certain volume from 4D dataset. For clarity, I added numbers to individual dots indicating slice numbers in figure below.

shouldn't the simulated data be constructed by fitting one polynomial function to the odd slices, and another one to the even slices?

You mean something like is plotted in the figure below? Then, however, individual pixel shift (blue dots) have to vary for odd and even slices. Otherwise, the fit would be linear (a.). I plotted two scenarios of variation for odd and even slices: constant variation (b.) or random variation (controlled by a factor changing in some interval) (c.).

[valosekj]

I generated some dummy 4D interleaved data using two polynomials (one fitting even slices and another one fitting odd slices) based on the previous comment. Voxel shift/motion varies randomly in interval from 0.5 to 3.5 for all examples below.

Example of 3rd order polynomials:

Example of 3rd order polynomials after moco_interleaved (in red are changes caused by moco):

Example of fitting 3rd order polynomials for one volume:

Example of 10th order polynomials:

Example of 10th order polynomials after moco_interleaved (in red are changes caused by moco):

Example of fitting 10th order polynomials for one volume:

Higher order of polynomials causes larger variation in S-I direction. So, I think that slowly-varying changes in the B0 field are better captured by polynomials of lower order (e.g. 3rd). @jcohenadad what do you think?

[jcohenadad]

@valosekj i'm having hard time understanding the overall workflow. To me, the synthetic data could be made by:

• create two synthetic segmentations (e.g. the one produced by dummy_segmentation()), by choosing random scalar for each term of the polynomial function (i.e. y=sum{a_i*x^i}, with a_i ∈ |-1, 1[)
• interleave these two segmentations (one for odd, the other for even)
• loop across time points.

So, i don't really see the need for "fitting" polynomial functions here.

[valosekj]

• create two synthetic segmentations (e.g. the one produced by dummy_segmentation()), by choosing random scalar for each term of the polynomial function (i.e. y=sum{a_i*x^i}, with a_i ∈ |-1, 1[)

The current implementation of dummy_segmentation function does not allow polynomial function:

https://github.com/neuropoly/spinalcordtoolbox/blob/cbc18e92806138bd612957ab71b6e9cd394fa00d/spinalcordtoolbox/testing/create_test_data.py#L140-L152

I implemented polynomial function only for interleaved option, but not globally for dummy_segmentation. However, I implemented it in bad way (discussion above).
So, I'll start with reimplementation from scratch. The first step will be update of dummy_segmentation function to allow polynomial function, right? In one of your previous comment you suggested to use code from dummy_centrerline. But I can see now, that dummy_centrerline switched from polynomial function to regularized curve meanwhile in this PR #2754.

Instead of this, I can probably use polyfit_1d function to produce polynomial function, right?: https://github.com/neuropoly/spinalcordtoolbox/blob/cbc18e92806138bd612957ab71b6e9cd394fa00d/spinalcordtoolbox/centerline/curve_fitting.py#L15

[jcohenadad]

The first step will be update of dummy_segmentation function to allow polynomial function, right?

👍

In one of your previous comment you suggested to use code from dummy_centrerline. But I can see now, that dummy_centrerline switched from polynomial function to regularized curve meanwhile in this PR #2754. Instead of this, I can probably use polyfit_1d function to produce polynomial function, right?

👍

[jcohenadad]

@valosekj are you still working on that PR?

[valosekj]

@valosekj are you still working on that PR?

@jcohenadad very sorry, when I worked on #2835 I somehow forgot to finish and merge this PR.
I can see now that two script modified by this PR (sct_dmri_moco.py and moco.py) have been moved meanwhile:

• scripts/sct_dmri_moco.py --> spinalcordtoolbox/spinalcordtoolbox/scripts/sct_dmri_moco.py

• spinalcordtoolbox/moco.py --> spinalcordtoolbox/spinalcordtoolbox/moco.py

I am not sure now, if I should try to rebase this PR or start from scratch with moved scripts.

[jcohenadad]

I am not sure now, if I should try to rebase this PR or start from scratch with moved scripts.

the scripts only moved, so you should be able to resolve conflicts easily. What I would do (i know that's not the proper git way), would be to

• copy the files you worked on in a temp folder outside of SCT dir,
• git merge master (to update your branch with master)
• copy back the files you worked on in the proper directory

@kousu @joshuacwnewton might have other/better advices

[valosekj]

What I would do (i know that's not the proper git way), would be to

• copy the files you worked on in a temp folder outside of SCT dir,
• git merge master (to update your branch with master)
• copy back the files you worked on in the proper directory

Thanks, I have tried it but I have had troubles with unresolved conflicts and unsuccessful rebase.
I have created locally new branch derived from current master and implemented interleaved option for both sct_dmri_moco.py and moco.py again. If it is not problem, I would close this PR and discard this branch and open new PR from my fresh new branch (which is local so far).

[valosekj]

Closed due to out-of-date with master and unresolved conflicts. Replaced by #3172