SDC-pepolar: topup implementation v2

[markushs]

Ref #37

Decided to do this as a new PR as the initial comments to #76 required implementing the approach somewhat differently.

pepolar.py is now modified to run fsl topup. pepolar.py returns a topup estimated fieldmap (in Hz) and a "magnitude" image (the unwarped epis from topup).

Next, the topup estimated fieldmap is converted to a displacement field using routines from the sdcflows.workflows.fmap module (fmap_wf + fmap2field_wf).

Finally, unwarping is performed using the sdcflow_workflows.unwarp module.

[pull-assistant]

Best reviewed: commit by commit

---

Optimal code review plan

     pepolar.py now contains topup workflow

     fix: fsl.Info().version() returns string

     connect pepolar wf with FIELDMAP path

Powered by Pull Assistant. Last update c5463e7 ... 44f1a3f. Read the comment docs.

[markushs]

Hi @oesteban. Is there still any interest in this approach? I'm asking as it was mentioned in a SDC-issue over at fmriprep. I had the impression that you were working on an alternative, and likely improved, implementation under dmriprep. If that's the case, should I just close the draft?

[oesteban]

Hi @markushs, thanks for bringing this up. The situation, indeed, is a bit messy.

We do have a functional alternative working for dMRIPrep -> https://github.com/nipreps/dmriprep/blob/a2cd0decb21e18697bec56942ec9c6dbda7e7e1b/dmriprep/workflows/fmap/base.py

While developing that, I realized of a few issues that are holding me back in upstreaming that code on to SDCFlows (and hence, make it available to fMRIPrep):

• The need for a more sophisticated specification from BIDS as to what was the intent of the researcher for fieldmapping, when she/he designed the acquisition protocol ([ENH] Allow encoding the fieldmapping intent of the protocol bids-standard/bids-specification#622). In the end, that would resolve (or alleviate) a longstanding issue for SDCFlows (RFC: New heuristics table #101).
• TOPUP from FSL 5 is quite unreliable and crashes with Segmentation Faults for a number of reasons (which are not undisclosed with the error). In dMRIPrep we basically patch TOPUP from FSL 6 into the overall FSL 5 installation. This solution is definitely not transparent (and probably unacceptable for fMRIPrep - opinions @effigies @mgxd ?). Therefore, any implementation of TOPUP for fMRIPrep entails upgrading FSL to v6. This is heavily problematic as Nipype itself does not cope with FSL 6 very smoothly.

I think both issues (BIDS specs and FSL 6) are the two main blockers for this. Both have workarounds, but I honestly think they are not transparent nor worth our time. But happy to be convinced on the contrary or learning of a more creative approach.

[oesteban]

We've made some progress in uniformizing the estimation API, and now we even have a base implementation for TOPUP:

sdcflows/sdcflows/models/pepolar.py

Lines 25 to 120 in 797eccb

	def init_topup_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"):
	"""
	Create the PEPOLAR field estimation workflow based on FSL's ``topup``.
	Workflow Graph
	.. workflow ::
	:graph2use: orig
	:simple_form: yes
	from sdcflows.models.pepolar import init_topup_wf
	wf = init_topup_wf()
	Parameters
	----------
	debug : :obj:`bool`
	Whether a fast configuration of topup (less accurate) should be applied.
	name : :obj:`str`
	Name for this workflow
	omp_nthreads : :obj:`int`
	Parallelize internal tasks across the number of CPUs given by this option.
	Inputs
	------
	in_data : :obj:`list` of :obj:`str`
	A list of EPI files that will be fed into TOPUP.
	metadata : :obj:`list` of :obj:`dict`
	A list of dictionaries containing the metadata corresponding to each file
	in ``in_data``.
	Outputs
	-------
	fmap : :obj:`str`
	The path of the estimated fieldmap.
	fmap_ref : :obj:`str`
	The path of an unwarped conversion of files in ``in_data``.
	"""
	from nipype.interfaces.fsl.epi import TOPUP
	from niworkflows.interfaces.nibabel import MergeSeries
	from sdcflows.interfaces.epi import GetReadoutTime
	workflow = Workflow(name=name)
	workflow.__desc__ = f"""\
	{_PEPOLAR_DESC} with `topup` @topup (FSL {TOPUP().version}).
	"""
	inputnode = pe.Node(
	niu.IdentityInterface(fields=["metadata", "in_data"]), name="inputnode"
	)
	outputnode = pe.Node(
	niu.IdentityInterface(
	fields=[
	"fmap_ref",
	"fmap",
	"coefficients",
	"jacobians",
	"xfms",
	"out_warps",
	]
	),
	name="outputnode",
	)
	concat_blips = pe.Node(MergeSeries(), name="concat_blips")
	readout_time = pe.MapNode(
	GetReadoutTime(),
	name="readout_time",
	iterfield=["metadata", "in_file"],
	run_without_submitting=True,
	)
	topup = pe.Node(
	TOPUP(
	config=_pkg_fname("sdcflows", f"data/flirtsch/b02b0{'_quick' * debug}.cnf")
	),
	name="topup",
	)
	# fmt: off
	workflow.connect([
	(inputnode, concat_blips, [("in_data", "in_files")]),
	(inputnode, readout_time, [("in_data", "in_file"),
	("metadata", "metadata")]),
	(inputnode, topup, [(("metadata", _pe2fsl), "encoding_direction")]),
	(readout_time, topup, [("readout_time", "readout_times")]),
	(concat_blips, topup, [("out_file", "in_file")]),
	(topup, outputnode, [("out_corrected", "fmap_ref"),
	("out_field", "fmap"),
	("out_fieldcoef", "coefficients"),
	("out_jacs", "jacobians"),
	("out_mats", "xfms"),
	("out_warps", "out_warps")]),
	])
	# fmt: on
	return workflow

I'm sure there are many edge cases we will need to account for, but the foundation is there. We even have a test running TOPUP: https://697-189292208-gh.circle-artifacts.com/0/tmp/tests/sdcflows/sub-HCP101006/figures/sub-HCP101006_desc-pepolar_fieldmap.svg

The idea is that all estimation strategies (namely, B0 fieldmaps acquired with Spiral Echo Imaging, phase-difference fieldmaps --both 2-echos and phasediff--, pepolar-AFNI, pepolar-TOPUP, end fieldmapless-SyN) produce not just a processed version of the fieldmap, also a field of B-Splines coefficients (i.e., the fieldcoeff output of TOPUP) -- please follow progress on this at #119.

This goes along with a more clear separation of estimation and correction. Estimation will generate the field and coefficients. These can be moved to the target image to be corrected, the field interpolated from the coefficients in the new space and finally, scaled by the Readout Time to generate the actual displacements field.

[oesteban]

Hi @markushs I think we should redirect these efforts towards the implementation we have on the dev/1.4.0 branch. Please send us all the feedback you have on that branch. Also, please let me know if you want to become a maintainer of SDCFlows, so that you can more flexibly contribute and we can also assign you reviews of PRs. We'd be excited that you joined forces with the project. WDYT?

In regards this particular PR, can I close it?

[markushs]

Hi @oesteban
Great that you seem to have found a better solution for this! I'll pay attention to the dev-branch from now on (no prob closing this PR). Concerning your question, I am flattered, but not sure if I can contribute with much. My python skills are getting better though (after a long and painful transition from matlab), so happy to help if there's anything I can do.

[oesteban]

Concerning your question, I am flattered, but not sure if I can contribute with much. My python skills are getting better though (after a long and painful transition from matlab), so happy to help if there's anything I can do.

You have already contributed with much - the fact that your two PRs on this TOPUP implementation haven't made it to master doesn't mean they don't have contributed to SDCFlows overall. Your initiative has stimulated the software go forward and has given ideas to others about how to proceed to make progress.

Looking into the future, although you might not speak Python fluently today I'm sure your feedback will be very valuable on others' PRs. And that's only the start.

Please let me know if you want to be included - as I said, we'd be honored :).

Closing this for the time being.