From 0ad307b66e5c97e07aa4ca0acf43abaeb4e3431e Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Wed, 20 Oct 2021 11:05:34 +0200
Subject: [PATCH] fix: ME | Revise SDC compute graph

Although I brought up the issue within #2530 and got some mildly opposed initial
opinion from @effigies
(https://github.com/nipreps/fmriprep/pull/2530#issuecomment-947147250),
I've come to think that we should provide the T2star workflow with
HMC'ed and SDC'ed echos (confirmation awaiting:
https://github.com/nipreps/fmriprep/pull/2608#pullrequestreview-784145171).

This comes as an issue related to #2606.

This PR proposes that individual echoes are HMC'ed (and SDC'ed when
fieldmaps are available) before calculating the optimal combination
and the T2star map.

I believe this may be beneficial when we combine SDC with modulation by
the Jacobian of the distortion (nipreps/sdcflows#238) because that might
help restore some (little, admittedly) of the dropout of the later
echoes in distorted regions. I would expect more voxels will be deemed
as acceptable for the exponential fitting for this reason.

It also conceptually simplifies a little, as both SE and ME paths look
more alike, and only the little hacks to generate ME iterables and the
``boldbuffer`` are now necessary. Otherwise, SDC would be inserted at a
much later stage for ME only.
---
 fmriprep/workflows/bold/base.py | 25 ++++++++++++-------------
 1 file changed, 12 insertions(+), 13 deletions(-)

diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index 684691c79..6ec411624 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -488,14 +488,15 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
                 (meepi_echos, bold_stc_wf, [("bold_file", "inputnode.bold_file")]),
             ])
             # fmt:on
-    elif not multiecho:  # STC is too short or False
+
+    # bypass STC from original BOLD in both SE and ME cases
+    elif not multiecho:  # SE and skip-STC
         # fmt:off
-        # bypass STC from original BOLD to the splitter through boldbuffer
         workflow.connect([
             (initial_boldref_wf, boldbuffer, [("outputnode.bold_file", "bold_file")]),
         ])
         # fmt:on
-    else:
+    else:  # ME and skip-STC
         # for meepi, iterate over all meepi echos to boldbuffer
         boldbuffer.iterables = ("bold_file", bold_file)
 
@@ -999,9 +1000,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
                 ]),
             ] if not multiecho
             else [
-                (bold_bold_trans_wf, join_echos, [
-                    ("outputnode.bold", "bold_files"),
-                ]),
                 (join_echos, final_boldref_wf, [("bold_files", "inputnode.bold_file")]),
                 # use reference image mask used by bold_bold_trans_wf
                 (bold_bold_trans_wf, bold_t2s_wf, [
@@ -1087,6 +1085,8 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
             ("outputnode.xforms", "inputnode.hmc_xforms")]),
         (initial_boldref_wf, sdc_report, [
             ("outputnode.ref_image", "before")]),
+        (bold_split, unwarp_wf, [
+            ("out_files", "inputnode.distorted")]),
         (unwarp_wf, final_boldref_wf, [
             ("outputnode.corrected_ref", "inputnode.bold_file"),
         ]),
@@ -1096,12 +1096,11 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
         (inputnode, ds_report_sdc, [("bold_file", "source_file")]),
         (sdc_report, ds_report_sdc, [("out_report", "in_file")]),
         # remaining workflow connections
-        (unwarp_wf, bold_std_trans_wf, [
+        (unwarp_wf, bold_t1_trans_wf, [
             # TEMPORARY: For the moment we can't use frame-wise fieldmaps
             (("outputnode.fieldwarp", _pop), "inputnode.fieldwarp"),
         ]),
-        (unwarp_wf, bold_final, [("outputnode.corrected", "bold")]),
-        (unwarp_wf, bold_t1_trans_wf, [
+        (unwarp_wf, bold_std_trans_wf, [
             # TEMPORARY: For the moment we can't use frame-wise fieldmaps
             (("outputnode.fieldwarp", _pop), "inputnode.fieldwarp"),
         ]),
@@ -1111,15 +1110,15 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
     if not multiecho:
         # fmt:off
         workflow.connect([
-            (bold_split, unwarp_wf, [
-                ("out_files", "inputnode.distorted")]),
+            (unwarp_wf, bold_final, [("outputnode.corrected", "bold")]),
         ])
         # fmt:on
     else:
         # fmt:off
         workflow.connect([
-            (split_opt_comb, unwarp_wf, [
-                ("out_files", "inputnode.distorted")])
+            (unwarp_wf, join_echos, [
+                ("outputnode.outputnode", "bold_files"),
+            ]),
         ])
         # fmt:on