@oesteban In order to build a DOCMA workflow, I'm trying to make a node that splits 4D files into 3D ones, and I'd prefer to have a function that flexibly allows any inputs so I don't need a separate node for each file I want to split.

def split_files(*, volume_idx, **kwargs):
    """
    Split an input 4D file and write out the set of 3D files.
    """
    import os.path as op
    from nilearn.image import index_img
    from nipype.utils.filemanip import split_filename

    out_files = []
    for in_file in kwargs.values():
        _, base, _ = split_filename(in_file)
        out_file = op.abspath(base + "_{0:05d}.nii.gz".format(volume_idx))
        img = index_img(in_file, volume_idx)
        img.to_filename(out_file)
        out_files.append(out_file)
    return out_files

With the above, I can do the following to make my splitting node:

buffernode = pe.Node(niu.IdentityInterface(
    fields=['file1_4d', 'file2_4d', 'file3_4d', 'volume_idx']),
    iterables=[('volume_idx', np.arange(num_trs, dtype=int))],
    name='buffernode')

split_4d_to_3d = pe.Node(
    interface=niu.Function(['file1_4d', 'file2_4d', 'file3_4d'],
                           ['files1_3d', 'files2_3d', 'files3_3d'],
                           split_files),
    name='split_4d_to_3d')

However, the function split_files requires that the kwargs retain their order, which is only implemented in Python 3.6+, so my question is: do you know of a good Python 3.5+ alternative or, if not, is it preferable to bump up requirements to Python 3.6+ or to have repetitive nodes where the function only splits one file and there's a separate node in the workflow for each file?

Maybe:

buffernode = pe.Node(niu.IdentityInterface(
    fields=['file1_4d', 'file2_4d', 'file3_4d', 'volume_idx']),
    iterables=[('volume_idx', np.arange(num_trs, dtype=int))],
    name='buffernode')

split_4d_to_3d = pe.MapNode(
    interface=niu.Function(['file', 'volume_idx'], ['file'], split_files),
    name='split_4d_to_3d', iterfield=['file', 'volume_idx'])

join_files = pe.JoinNode(...)

@oesteban Okay I have a workflow that successfully runs, although the unwarping results are nonsensical (see below). Should I open a draft PR so you can review the code directly or would you prefer to continue discussion here?

First volume's field map

First volume unwarped

I'm trying out ROMEO for the unwrapping, and possibly even B0-map generation, so hopefully I'll be able to move forward on this soon.

@vanandrew has developed a very cool dynamic distortion method for complex-valued, multi-echo BOLD data called MEDIC (https://github.com/vanandrew/medic_analysis). There's a new preprint about it available here, and it's integrated into the Dosenbach lab's preprocessing pipeline (https://github.com/DosenbachGreene/processing_pipeline). I think it can be safely considered a superior alternative to DOCMA.

EDIT: One additional consideration is that dynamic distortion correction would probably be best applied before HMC, unlike static distortion correction.

Since there's the "gathering feedback" tag here, I just wanted throw my two cents in and say that the framewise SDC from Andrew Van (@vanandrew) seems like it would be a great addition to FMRIprep, but I am not sure how it interacts with the Echo Time Difference check in the BIDS validator. In our initial tests, we couldn't get this to run since our echoes are about 17 ms apart and the BIDS validator looks like it's expecting a max of 10 ms. More here: https://neurostars.org/t/multi-echo-fmaps-echo-1-and-2-difference-unreasonable/28846

Thanks!
David

The BIDS validator check can be relaxed if there are false positives. The idea was just to warn people they almost certainly made a mistake during conversion or curation.

FWIW I don't think SDCFlows will work if you just take the first two echoes from one of your volumes and rename them as field map files (which is when the BIDS validator would notice your echo time difference, I believe). The phase data are probably highly wrapped, which I believe is what ROMEO is particularly good at handling, but ROMEO isn't what SDCFlows uses for unwrapping. I could definitely be wrong though.

Yes, there's definitely a ton of wrapping with the phase images. I'm not yet 100% sure what MEDIC is doing under the hood, but the output we got was quite reasonable at least in terms of how the images looked.

We had four echoes going into MEDIC, but it doesn't output a corresponding .json file with the metadata (and I was assuming FMRIprep was using the .json file to infer parameters about the scan, like the echo time difference). In the case of a multi-echo sequence with four echoes where MEDIC produces a great fmap image for each time point, I'm just not sure how FMRIprep should handle the fmap. If we had a proper .json file for any of these fmaps (ideally the first one, though), it seems like FMRIprep could quickly and easily have task- and run-specific fmaps that could be applied to each bold image. Of course, I realize that the major advance with Andrew's paper is doing this in a framewise fashion, but I know that will take more time.