[WIP] Eddy Motion Correction (EMC) #62
Conversation
|
Best reviewed: commit by commit
Optimal code review plan (4 warnings)
|
|
Hello @dPys, Thank you for updating!
To test for issues locally, Comment last updated at 2020-03-24 17:24:31 UTC |
Codecov Report
@@ Coverage Diff @@
## master #62 +/- ##
==========================================
- Coverage 52.37% 42.17% -10.2%
==========================================
Files 16 21 +5
Lines 907 1648 +741
Branches 114 172 +58
==========================================
+ Hits 475 695 +220
- Misses 431 952 +521
Partials 1 1
Continue to review full report at Codecov.
|
dPys
force-pushed
the
HMC
branch
2 times, most recently
from
dab9f38 to
9045994
Compare
| "transform_parameters": [ [ 0.2 ], [ 0.15 ] ], | ||
| "convergence_threshold": [ 1e-06, 1e-06 ], | ||
| "convergence_window_size": [ 20, 20 ], | ||
| "metric": [ "Mattes", "Mattes" ], |
There was a problem hiding this comment.
Is this choice based on speed?
Wouldn't it make more sense some cross-correlation based measure? MI is just going to exacerbate the problems you mentioned about the different content of each direction (pinging @arokem and @mattcieslak).
There was a problem hiding this comment.
I agree @oesteban -- cross-corr seems like it would be preferable here. Curious what @mattcieslak thinks?
There was a problem hiding this comment.
Remember, the fixed (model-estimated) images here will have the same contrast as the moving images so most metrics would be ok here. MI is much faster than CC, so any benefit of CC would be outweighed by its extra computation time. Also, this is only the coarse iteration.
There was a problem hiding this comment.
Yup, but MI has much flatter gradient - I'm going to make up a number here, but 1 CC iteration is as effective as 10 MI iterations.
dmriprep/interfaces/bids.py
Outdated
|
|
||
|
|
||
| class DerivativesDataSink(SimpleInterface): | ||
| """ | ||
| Saves the `in_file` into a BIDS-Derivatives folder provided | ||
| by `base_directory`, given the input reference `source_file`. | ||
| >>> from pathlib import Path | ||
| >>> import tempfile | ||
| >>> from dmriprep.utils.bids import collect_data | ||
| >>> tmpdir = Path(tempfile.mkdtemp()) | ||
| >>> tmpfile = tmpdir / 'a_temp_file.nii.gz' | ||
| >>> tmpfile.open('w').close() # "touch" the file | ||
| >>> dsink = DerivativesDataSink(base_directory=str(tmpdir)) | ||
| >>> dsink.inputs.in_file = str(tmpfile) | ||
| >>> dsink.inputs.source_file = collect_data('ds114', '01')[0]['t1w'][0] | ||
| >>> dsink.inputs.keep_dtype = True | ||
| >>> dsink.inputs.suffix = 'target-mni' | ||
| >>> res = dsink.run() | ||
| >>> res.outputs.out_file # doctest: +ELLIPSIS | ||
| '.../dmriprep/sub-01/ses-retest/anat/sub-01_ses-retest_target-mni_T1w.nii.gz' | ||
| >>> bids_dir = tmpdir / 'bidsroot' / 'sub-02' / 'ses-noanat' / 'func' | ||
| >>> bids_dir.mkdir(parents=True, exist_ok=True) | ||
| >>> tricky_source = bids_dir / 'sub-02_ses-noanat_task-rest_run-01_bold.nii.gz' | ||
| >>> tricky_source.open('w').close() | ||
| >>> dsink = DerivativesDataSink(base_directory=str(tmpdir)) | ||
| >>> dsink.inputs.in_file = str(tmpfile) | ||
| >>> dsink.inputs.source_file = str(tricky_source) | ||
| >>> dsink.inputs.keep_dtype = True | ||
| >>> dsink.inputs.desc = 'preproc' | ||
| >>> res = dsink.run() | ||
| >>> res.outputs.out_file # doctest: +ELLIPSIS | ||
| '.../dmriprep/sub-02/ses-noanat/func/sub-02_ses-noanat_task-rest_run-01_\ | ||
| desc-preproc_bold.nii.gz' | ||
| """ | ||
| input_spec = DerivativesDataSinkInputSpec | ||
| output_spec = DerivativesDataSinkOutputSpec | ||
| out_path_base = "dmriprep" | ||
| _always_run = True | ||
|
|
||
| def __init__(self, out_path_base=None, **inputs): | ||
| super(DerivativesDataSink, self).__init__(**inputs) | ||
| self._results['out_file'] = [] | ||
| if out_path_base: | ||
| self.out_path_base = out_path_base | ||
|
|
||
| def _run_interface(self, runtime): | ||
| src_fname, _ = _splitext(self.inputs.source_file) | ||
| src_fname, dtype = src_fname.rsplit('_', 1) | ||
| _, ext = _splitext(self.inputs.in_file[0]) | ||
| if self.inputs.compress is True and not ext.endswith('.gz'): | ||
| ext += '.gz' | ||
| elif self.inputs.compress is False and ext.endswith('.gz'): | ||
| ext = ext[:-3] | ||
|
|
||
| m = BIDS_NAME.search(src_fname) | ||
|
|
||
| mod = os.path.basename(os.path.dirname(self.inputs.source_file)) | ||
|
|
||
| base_directory = runtime.cwd | ||
| if isdefined(self.inputs.base_directory): | ||
| base_directory = str(self.inputs.base_directory) | ||
|
|
||
| out_path = '{}/{subject_id}'.format(self.out_path_base, **m.groupdict()) | ||
| if m.groupdict().get('session_id') is not None: | ||
| out_path += '/{session_id}'.format(**m.groupdict()) | ||
| out_path += '/{}'.format(mod) | ||
|
|
||
| out_path = os.path.join(base_directory, out_path) | ||
|
|
||
| os.makedirs(out_path, exist_ok=True) | ||
|
|
||
| if isdefined(self.inputs.prefix): | ||
| base_fname = os.path.join(out_path, self.inputs.prefix) | ||
| else: | ||
| base_fname = os.path.join(out_path, src_fname) | ||
|
|
||
| formatstr = '{bname}{space}{desc}{suffix}{dtype}{ext}' | ||
| if len(self.inputs.in_file) > 1 and not isdefined(self.inputs.extra_values): | ||
| formatstr = '{bname}{space}{desc}{suffix}{i:04d}{dtype}{ext}' | ||
|
|
||
| space = '_space-{}'.format(self.inputs.space) if self.inputs.space else '' | ||
| desc = '_desc-{}'.format(self.inputs.desc) if self.inputs.desc else '' | ||
| suffix = '_{}'.format(self.inputs.suffix) if self.inputs.suffix else '' | ||
| dtype = '' if not self.inputs.keep_dtype else ('_%s' % dtype) | ||
|
|
||
| self._results['compression'] = [] | ||
| for i, fname in enumerate(self.inputs.in_file): | ||
| out_file = formatstr.format( | ||
| bname=base_fname, | ||
| space=space, | ||
| desc=desc, | ||
| suffix=suffix, | ||
| i=i, | ||
| dtype=dtype, | ||
| ext=ext) | ||
| if isdefined(self.inputs.extra_values): | ||
| out_file = out_file.format(extra_value=self.inputs.extra_values[i]) | ||
| self._results['out_file'].append(out_file) | ||
| self._results['compression'].append(_copy_any(fname, out_file)) | ||
| return runtime |
There was a problem hiding this comment.
We already have a DerivativesDataSink:
dmriprep/dmriprep/interfaces/__init__.py
Lines 13 to 16 in 0a609b4
Unless there is some new functionality here, this can go away.
There was a problem hiding this comment.
yep! (didn't see that it already lived in init.py initially.
dmriprep/interfaces/images.py
Outdated
| class _RescaleB0InputSpec(BaseInterfaceInputSpec): | ||
| in_file = File(exists=True, mandatory=True, desc='b0s file') | ||
| mask_file = File(exists=True, mandatory=True, desc='mask file') | ||
|
|
||
|
|
||
| class _RescaleB0OutputSpec(TraitedSpec): | ||
| out_ref = File(exists=True, desc='One average b0 file') | ||
| out_ref = File(exists=True, desc='One median b0 file') |
There was a problem hiding this comment.
Please avoid changes out of scope. This one should've been picked up during the review of #25 - that said, I think it is better "average" here, as it needs to be informative.
dmriprep/interfaces/images.py
Outdated
| class ImageMath(CommandLine): | ||
| input_spec = ImageMathInputSpec | ||
| output_spec = ImageMathOutputSpec | ||
| _cmd = 'ImageMath' | ||
|
|
||
| def _gen_filename(self, name): | ||
| if name == 'out_file': | ||
| output = self.inputs.out_file | ||
| if not isdefined(output): | ||
| _, fname, ext = split_filename(self.inputs.in_file) | ||
| output = fname + "_" + self.inputs.operation + ext | ||
| return output | ||
| return None | ||
|
|
||
| def _list_outputs(self): | ||
| outputs = self.output_spec().get() | ||
| outputs['out_file'] = os.path.abspath(self._gen_filename('out_file')) | ||
| return outputs |
There was a problem hiding this comment.
Please use ImageMath from niworkflows (and PR against niworkflows if some functionality is missing)
There was a problem hiding this comment.
So the reason for redefining it outside of niworkflows is that we need it to be able to take an explicit secondary argument and file. In other words, as defined here, it is actually quite different. But we could monkey-patch it?
There was a problem hiding this comment.
Yes, we can import the niworkflows version and then add to it (similar to the DerivativesDataSink and BIDSDataGrabberOutputSpec examples here: https://github.com/nipreps/dmriprep/blob/master/dmriprep/interfaces/__init__.py#L13-L20)
dmriprep/interfaces/vectors.py
Outdated
| class ReorientVectors(SimpleInterface): | ||
| """ | ||
| Reorient Vectors | ||
| Example | ||
| ------- | ||
| >>> os.chdir(tmpdir) | ||
| >>> oldrasb = str(data_dir / 'dwi.tsv') | ||
| >>> oldrasb_mat = np.loadtxt(str(data_dir / 'dwi.tsv'), skiprows=1) | ||
| >>> # The simple case: all affines are identity | ||
| >>> affine_list = np.zeros((len(oldrasb_mat[:, 3][oldrasb_mat[:, 3] != 0]), 4, 4)) | ||
| >>> for i in range(4): | ||
| >>> affine_list[:, i, i] = 1 | ||
| >>> reor_vecs = ReorientVectors() | ||
| >>> reor_vecs = ReorientVectors() | ||
| >>> reor_vecs.inputs.affines = affine_list | ||
| >>> reor_vecs.inputs.in_rasb = oldrasb | ||
| >>> res = reor_vecs.run() | ||
| >>> out_rasb = res.outputs.out_rasb | ||
| >>> out_rasb_mat = np.loadtxt(out_rasb, skiprows=1) | ||
| >>> npt.assert_equal(oldrasb_mat, out_rasb_mat) | ||
| True | ||
| """ | ||
|
|
||
| input_spec = _ReorientVectorsInputSpec | ||
| output_spec = _ReorientVectorsOutputSpec | ||
|
|
||
| def _run_interface(self, runtime): | ||
| from nipype.utils.filemanip import fname_presuffix | ||
| reor_table = reorient_vecs_from_ras_b( | ||
| rasb_file=self.inputs.rasb_file, | ||
| affines=self.inputs.affines, | ||
| b0_threshold=self.inputs.b0_threshold, | ||
| ) | ||
|
|
||
| cwd = Path(runtime.cwd).absolute() | ||
| reor_rasb_file = fname_presuffix( | ||
| self.inputs.rasb_file, use_ext=False, suffix='_reoriented.tsv', | ||
| newpath=str(cwd)) | ||
| np.savetxt(str(reor_rasb_file), reor_table, | ||
| delimiter='\t', header='\t'.join('RASB'), | ||
| fmt=['%.8f'] * 3 + ['%g']) | ||
|
|
||
| self._results['out_rasb'] = reor_rasb_file | ||
| return runtime |
There was a problem hiding this comment.
Yes, but with added functionality. All future changes for this interface will be redirected to the #49 PR
| def splitext(fname): | ||
| """Splits filename and extension (.gz safe) | ||
| >>> splitext('some/file.nii.gz') | ||
| ('file', '.nii.gz') | ||
| >>> splitext('some/other/file.nii') | ||
| ('file', '.nii') | ||
| >>> splitext('otherext.tar.gz') | ||
| ('otherext', '.tar.gz') | ||
| >>> splitext('text.txt') | ||
| ('text', '.txt') | ||
| """ | ||
| from pathlib import Path | ||
| basename = str(Path(fname).name) | ||
| stem = Path(basename.rstrip('.gz')).stem | ||
| return stem, basename[len(stem):] | ||
|
|
||
|
|
||
| def _copy_any(src, dst): | ||
| import os | ||
| import gzip | ||
| from shutil import copyfileobj | ||
| from nipype.utils.filemanip import copyfile | ||
|
|
||
| src_isgz = src.endswith('.gz') | ||
| dst_isgz = dst.endswith('.gz') | ||
| if not src_isgz and not dst_isgz: | ||
| copyfile(src, dst, copy=True, use_hardlink=True) | ||
| return False # Make sure we do not reuse the hardlink later | ||
|
|
||
| # Unlink target (should not exist) | ||
| if os.path.exists(dst): | ||
| os.unlink(dst) | ||
|
|
||
| src_open = gzip.open if src_isgz else open | ||
| with src_open(src, 'rb') as f_in: | ||
| with open(dst, 'wb') as f_out: | ||
| if dst_isgz: | ||
| # Remove FNAME header from gzip (poldracklab/fmriprep#1480) | ||
| gz_out = gzip.GzipFile('', 'wb', 9, f_out, 0.) | ||
| copyfileobj(f_in, gz_out) | ||
| gz_out.close() | ||
| else: | ||
| copyfileobj(f_in, f_out) | ||
|
|
||
| return True |
There was a problem hiding this comment.
I think this comes from nipype, why not just import them?
There was a problem hiding this comment.
Didn't see these functions in nipype, but maybe I missed them somewhere?
| @@ -0,0 +1,113 @@ | |||
| import numpy as np | |||
There was a problem hiding this comment.
This new file makes a lot of sense, especially if we can come up with unit tests for all these helper functions.
Could you send a separate PR to split the currently existing functions (extract_b0 and rescale_b0) out to this file so they don't show up in the diff?
dmriprep/utils/vectors.py
Outdated
| @@ -375,3 +375,109 @@ def bvecs2ras(affine, bvecs, norm=True, bvec_norm_epsilon=0.2): | |||
| rotated_bvecs[~b0s] /= norms_bvecs[~b0s, np.newaxis] | |||
| rotated_bvecs[b0s] = np.zeros(3) | |||
| return rotated_bvecs | |||
|
|
|||
|
|
|||
| def reorient_vecs_from_ras_b(rasb_file, affines, b0_threshold=B0_THRESHOLD): | |||
There was a problem hiding this comment.
I would try to keep the b0 threshold encoded in the vector representation object, meaning, no need to keep dragging it all around. Please remove.
There was a problem hiding this comment.
Actually, I would incorporate the reorientation into the vectors object. Will comment on this at #49. We definitely should push that one first.
|
Thank you @oesteban for all of your helpful feedback. I will be out of town next week so won't be able to work a ton on this after tomorrow. But will do my best to keep the dev of EMC moving... In the meantime, @arokem -- anything that we can think of to speed up the way that sfm predicts signal (i.e. from a given vector coordinate) would be helpful. If it can't be the ElasticNet itself, then we might have to get creative to make it viable with LOO here. Please let me know if there's anything that can be done. More soon. |
|
Yes. That's a good idea and my hunch is that it would work well: Ridge SFM!
Not too hard to implement either as the DIPY SFM object can take a solver
as an input.
…On Fri, Jan 17, 2020 at 4:40 PM Matt Cieslak ***@***.***> wrote:
@arokem <https://github.com/arokem> @dPys <https://github.com/dPys> one
thing that made shoreline reasonably fast with 3dSHORE is that I used the
L2 (with LShore and NShore) fit instead of L1 or elasticnet. This reduces
the problem to a couple dot products and still resulted in very good signal
prediction. Would this work with a sfm?
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#62?email_source=notifications&email_token=AAA46NUS7S252VRARTY2ERTQ6JFXBA5CNFSM4KHNRSUKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEJJLQ2Q#issuecomment-575846506>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AAA46NVYCQLNKKKKQFIKZ23Q6JFXBANCNFSM4KHNRSUA>
.
|
Rename to EMC (Eddy Motion Correction) and fix several workflow bugs
dPys
force-pushed
the
HMC
branch
2 times, most recently
from
cf5d7c2 to
8b9c73e
Compare
|
@arokem @oesteban @mattcieslak -- The EMC WIP should now run from start to finish. I've also gone ahead and replaced every ANTs reg call with a dipy / pure-python equivalent. These new reg routines are based on @arokem 's api.py file, but have been modified to read and write to files from disk. The routines have also been more or less tuned to the original coarse/precise rigid/affine reg parameters used with the ANTs calls. Another note-- With these further upgrades, I'm also seeing a dramatic speedup for EMC as a whole. To improve efficiency, I initialized the signal prediction step using the tensor model on the first iteration, which allows us to get a really good coarse transformation to start. The subsequent iteration(s), used for precise transforms, are then performed with the more computationally expensive sfm model. At least from visual inspection, the impact of SFM for predicting signal with single-shell data is pretty excellent. A note: sfm is now incorporated with a ridge/L2 solver, which made prediction much faster. To get a sense of EMC's effectiveness at this point, here is a before-and-after:
If anyone wish to try running this on your own from my fork, here's how you could do it (note that it assumes a rough median/mean b0 template has already been created and can be used as input): Now, to actually PR this beast, I propose the following six separate PR's (@oesteban -- please confirm):
|
Clean up motion plotting
|
Taking up from where @dPys left off:
I think we have pretty much all we need already merged. Will look back at this if something is missing.
@arokem is working on this.
I'll extract the code from this PR and separate on to a new one.
I think this is already done in part of in full by @josephmje. As for 1, will come back to this PR if anything is needed.
Let's leave this out for now, this is being worked out on a different repo.
Will extract this after 3 is merged. |
|
Okay, pushed Derek's branch (and his commits) to https://github.com/nipreps/dmriprep/tree/feat/wip/hmc-dpys, just to be sure we don't lose any lines. |
|
Closing for clarity. |
*Per @arokem 's request to PR the WIP, even though tests + docs are not yet included.
*After an initial glance from folks, plan is to close this WIP PR and eventually separate it into several PR's, perhaps per workflow (and associated interfaces/funcs) @oesteban ?
*About 50% of the code is adapted from fmriprep utilities and qsiprep's implementation of 3dShore, & around 50% is new.
*Supports 3dShore EMC for multishell data, but now implements it exclusively in Dipy as opposed to BrainSuite.
*More importantly, can implement Tensor and SparseFascicle model for LOO prediction on single-shell data!
*Registration routines currently use ANTs, but could probably be adapted to use to Dipy instead. In that case, they would most likely require careful retuning.
*Pictures to come, but initial visual checks seem to indicate at least decent quality performance (perhaps even reasonably good performance) as a starting benchmark.