[ENH] move biascorrect so it runs on resampled data by default (#527)

.circleci/DRBUDDI_SHORELine_epi.sh

@@ -33,7 +33,7 @@ ${QSIPREP_CMD} \
 	--sloppy \
 	--dwi-only \
 	--denoise-method none \
-	--dwi-no-biascorr \
+	--b1-biascorrect-stage none \
 	--output-space T1w \
 	--pepolar-method DRBUDDI \
 	--hmc-model none \

.circleci/DRBUDDI_eddy_rpe_series.sh

@@ -34,7 +34,7 @@ ${QSIPREP_CMD} \
 	--sloppy \
 	--dwi-only \
 	--denoise-method none \
-	--dwi-no-biascorr \
+	--b1_biascorrect_stage none \
 	--output-space T1w \
 	--pepolar-method DRBUDDI \
 	--eddy_config ${EDDY_CFG} \

.circleci/DSCSDSI.sh

@@ -19,7 +19,7 @@ Inputs:
  - DSCSDSI BIDS data (data/DSCSDSI_nofmap)
 
 DOC
-set +e 
+set +e
 # Setup environment and get data
 source ./get_data.sh
 TESTDIR=${PWD}
@@ -42,7 +42,7 @@ ${QSIPREP_CMD} \
    --sloppy --write-graph --use-syn-sdc \
    --force-syn \
    --output-space T1w \
-   --dwi-no-biascorr \
+   --b1_biascorrect_stage none \
    --hmc_model 3dSHORE \
    --hmc-transform Rigid \
    --shoreline_iters 1 \

.circleci/DSDTI_TOPUP.sh

@@ -41,6 +41,7 @@ ${QSIPREP_CMD} \
 	--sloppy \
 	--unringing-method mrdegibbs \
 	--output-space T1w \
+	--b1-biascorrect-stage legacy \
 	--recon-spec dsi_studio_gqi \
 	--eddy_config ${EDDY_CFG} \
 	--output-resolution 5 \

.circleci/DSDTI_nofmap.sh

@@ -43,6 +43,7 @@ ${QSIPREP_CMD} \
      --eddy-config ${EDDY_CFG} \
      --denoise-method none \
      --unringing-method rpg \
+     --b1-biascorrect-stage none \
      --sloppy \
      --output-space T1w \
      --output-resolution 5 \

.circleci/DSDTI_synfmap.sh

@@ -44,6 +44,7 @@ ${QSIPREP_CMD} \
      --eddy-config ${EDDY_CFG} \
      --sloppy \
 	--force-syn \
+	--b1-biascorrect-stage final \
      --denoise-method none \
 	--output-space T1w \
 	--output-resolution 5 \

.circleci/IntramodalTemplate.sh

@@ -40,11 +40,11 @@ QSIPREP_CMD=$(run_qsiprep_cmd ${BIDS_INPUT_DIR} ${OUTPUT_DIR})
 ${QSIPREP_CMD} \
 	 -w ${TEMPDIR} \
 	 --sloppy \
-	 --denoise-method none \
+	 --b1_biascorrect_stage none \
 	 --output-space T1w \
 	 --hmc_model none \
 	 --b0-motion-corr-to first \
 	 --output-resolution 5 \
 	 --intramodal-template-transform BSplineSyN \
 	 --intramodal-template-iters 2 \
-	 -vv 
+	 -vv

docs/preprocessing.rst

@@ -506,7 +506,7 @@ DWI preprocessing
                               dwi_denoise_window=5,
                               denoise_method='dwidenoise',
                               unringing_method='mrdegibbs',
-                              dwi_no_biascorr=False,
+                              b1_biascorr_stage='final',
                               no_b0_harmonization=False,
                               denoise_before_combining=True,
                               template='MNI152NLin2009cAsym',

qsiprep/cli/run.py

@@ -274,7 +274,16 @@ def get_parser():
     g_conf.add_argument(
         '--dwi-no-biascorr', '--dwi_no_biascorr',
         action='store_true',
-        help='skip b0-based dwi spatial bias correction')
+        help='DEPRECATED: see --b1-biascorr-stage')
+    g_conf.add_argument(
+        "--b1-biascorrect-stage", "--b1_biascorrect_stage",
+        action="store",
+        choices=["final", "none", "legacy"],
+        default="final",
+        help="Which stage to apply B1 bias correction. The default 'final' will "
+             "apply it after all the data has been resampled to its final space. "
+             "'none' will skip B1 bias correction and 'legacy' will behave consistent "
+             "with qsiprep < 0.17.")
     g_conf.add_argument(
         '--no-b0-harmonization', '--no_b0_harmonization',
         action='store_true',
@@ -969,7 +978,7 @@ def build_qsiprep_workflow(opts, retval):
         pepolar_method=opts.pepolar_method,
         dwi_denoise_window=opts.dwi_denoise_window,
         unringing_method=opts.unringing_method,
-        dwi_no_biascorr=opts.dwi_no_biascorr,
+        b1_biascorrect_stage=opts.b1_biascorrect_stage,
         no_b0_harmonization=opts.no_b0_harmonization,
         denoise_before_combining=not opts.denoise_after_combining,
         write_local_bvecs=opts.write_local_bvecs,

qsiprep/interfaces/dwi_merge.py

@@ -12,6 +12,7 @@
 from .fmap import get_distortion_grouping
 from ..workflows.dwi.util import _get_concatenated_bids_name
 LOGGER = logging.getLogger('nipype.workflow')
+MAX_COMBINED_SCANS = 100
 
 
 class MergeDWIsInputSpec(BaseInterfaceInputSpec):
@@ -214,6 +215,108 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _SplitResampledDWIsInputSpec(BaseInterfaceInputSpec):
+    dwi_file = File(exists=True, mandatory=True)
+    bval_file = File(exists=True, mandatory=True)
+    bvec_file = File(exists=True, mandatory=True)
+    confounds = File(exists=True, mandatory=True)
+    n_images = traits.Int(1)
+
+
+class _SplitResampledDWIsOutputSpec(TraitedSpec):
+    pass
+
+
+# Add slots for the possibly
+for subscan in np.arange(MAX_COMBINED_SCANS) + 1:
+    _SplitResampledDWIsOutputSpec.add_class_trait(
+        "dwi_file_%d" % subscan, File(exists=True))
+    _SplitResampledDWIsOutputSpec.add_class_trait(
+        "bval_file_%d" % subscan, File(exists=True))
+    _SplitResampledDWIsOutputSpec.add_class_trait(
+        "bvec_file_%d" % subscan, File(exists=True))
+    _SplitResampledDWIsOutputSpec.add_class_trait(
+        "source_file_%d" % subscan, traits.Str())
+
+
+class SplitResampledDWIs(SimpleInterface):
+    input_spec = _SplitResampledDWIsInputSpec
+    output_spec = _SplitResampledDWIsOutputSpec
+
+    def _run_interface(self, runtime):
+        # Load the confounds
+        confounds_df = pd.read_csv(self.inputs.confounds, sep="\t")
+        original_files = confounds_df['original_file'].unique().tolist()
+        if not len(original_files) == self.inputs.n_images:
+            raise Exception(
+                "Found %d files in confounds file, but expected %d" % (
+                    len(original_files), self.inputs.n_images))
+        resampled_img = load_img(self.inputs.dwi_file)
+        for file_num, original_file in enumerate(original_files, start=1):
+            image_indices = np.flatnonzero(
+                (confounds_df['original_file'] == original_file).to_numpy())
+            dwi_subfile = fname_presuffix(
+                original_file, prefix="resampled_", suffix=".nii.gz",
+                use_ext=False, newpath=runtime.cwd)
+            bval_subfile = dwi_subfile.replace(".nii.gz", ".bval")
+            bvec_subfile = dwi_subfile.replace(".nii.gz", ".bvec")
+            index_img(resampled_img, image_indices).to_filename(dwi_subfile)
+            subset_bvals(self.inputs.bval_file, image_indices, bval_subfile)
+            subset_bvecs(self.inputs.bvec_file, image_indices, bvec_subfile)
+
+            self._results["dwi_file_%d" % file_num] = dwi_subfile
+            self._results["bval_file_%d" % file_num] = bval_subfile
+            self._results["bvec_file_%d" % file_num] = bvec_subfile
+            self._results["source_file_%d" % file_num] = original_file
+        return runtime
+
+
+class _MergeFinalConfoundsInputSpec(BaseInterfaceInputSpec):
+    confounds = File(exists=True, mandatory=True)
+    bias_correction_confounds = InputMultiObject(File(exists=True), mandatory=False)
+    patch2self_correction_confounds = File(exists=True, mandatory=False)
+
+
+class _MergeFinalConfoundsOutputSpec(TraitedSpec):
+    confounds = File(exists=True)
+
+
+class MergeFinalConfounds(SimpleInterface):
+    input_spec = _MergeFinalConfoundsInputSpec
+    output_spec = _MergeFinalConfoundsOutputSpec
+
+    def _run_interface(self, runtime):
+
+        to_concat_horizontally = []
+        # New confounds from bias correction
+        if isdefined(self.inputs.bias_correction_confounds):
+            # There may be multuple files that need to be vertically stacked
+            biascorrection_df = pd.concat(
+                [pd.read_csv(bc_csv) for bc_csv in self.inputs.bias_correction_confounds],
+                axis=0, ignore_index=True)
+            to_concat_horizontally.append(biascorrection_df)
+        # New confounds from patch2self
+        if isdefined(self.inputs.patch2self_correction_confounds):
+            to_concat_horizontally.append(
+                pd.read_csv(self.inputs.patch2self_correction_confounds))
+
+        # If we have new ones, append the columns, prefixed by "final_"
+        if to_concat_horizontally:
+            new_confounds_file = fname_presuffix(
+                self.inputs.confounds, newpath=runtime.cwd, prefix="final_")
+            original_confounds = pd.read_csv(self.inputs.confounds, sep="\t")
+            extra_confounds = pd.concat(to_concat_horizontally, axis=1)
+            extra_confounds.columns = [
+                "final_" + col for col in extra_confounds.columns.tolist()]
+            final_confounds = pd.concat([original_confounds, extra_confounds], axis=1)
+            final_confounds.to_csv(new_confounds_file, sep="\t", index=False)
+            self._results['confounds'] = new_confounds_file
+        else:
+            self._results['confounds'] = self.inputs.confounds
+
+        return runtime
+
+
 def find_image_pairs(original_bvecs, bvals, assignments):
     assignments = np.array(assignments)
     group1_mask = assignments == 1
@@ -413,6 +516,18 @@ def _run_interface(self, runtime):
         return runtime
 
 
+def subset_bvals(bval_file, indices, out_bval_file):
+    original_bvals = np.loadtxt(bval_file)
+    bval_subset = original_bvals[indices]
+    np.savetxt(out_bval_file, bval_subset, fmt=str("%i"))
+
+
+def subset_bvecs(bvec_file, indices, out_bvec_file):
+    original_bvecs = np.loadtxt(bvec_file)
+    bvec_subset = original_bvecs[:, indices]
+    np.savetxt(out_bvec_file, bvec_subset, fmt=str("%.8f"))
+
+
 def combined_bval_array(bval_files):
     collected_vals = []
     for bval_file in bval_files:

qsiprep/interfaces/reports.py

@@ -498,10 +498,13 @@ def topup_selection_to_report(selected_indices, original_files, spec_lookup,
 
 
 class _SeriesQCInputSpec(BaseInterfaceInputSpec):
-    pre_qc = File(exists=True, desc='qc file from the raw data')
+    pre_qc = File(exists=True, desc='qc file from the raw data', mandatory=True)
     t1_qc = File(exists=True, desc='qc file from preprocessed image in t1 space')
-    mni_qc = File(exists=True, desc='qc file from preprocessed image in template space')
-    confounds_file = File(exists=True, desc='confounds file')
+    t1_qc_postproc = File(exists=True, desc='qc file from preprocessed image in template space')
+    confounds_file = File(exists=True, desc='confounds file', mandatory=True)
+    t1_mask_file = File(exists=True, desc='brain mask in t1 space')
+    t1_cnr_file = File(exists=True, desc='CNR file in t1 space')
+    t1_b0_series = File(exists=True, desc='time series of b=0 images')
     t1_dice_score = traits.Float()
     mni_dice_score = traits.Float()
     output_file_name = traits.File()
@@ -519,8 +522,8 @@ def _run_interface(self, runtime):
         image_qc = _load_qc_file(self.inputs.pre_qc, prefix="raw_")
         if isdefined(self.inputs.t1_qc):
             image_qc.update(_load_qc_file(self.inputs.t1_qc, prefix="t1_"))
-        if isdefined(self.inputs.mni_qc):
-            image_qc.update(_load_qc_file(self.inputs.mni_qc, prefix="mni_"))
+        if isdefined(self.inputs.t1_qc_postproc):
+            image_qc.update(_load_qc_file(self.inputs.t1_qc_postproc, prefix="t1post_"))
         motion_summary = calculate_motion_summary(self.inputs.confounds_file)
         image_qc.update(motion_summary)
 
@@ -530,6 +533,14 @@ def _run_interface(self, runtime):
         if isdefined(self.inputs.mni_dice_score):
             image_qc['mni_dice_distance'] = [self.inputs.mni_dice_score]
 
+        if isdefined(self.inputs.t1_mask_file):
+            if isdefined(self.inputs.t1_cnr_file):
+                # Add a function here to get CNR
+                pass
+            if isdefined(self.inputs.t1_b0_series):
+                # Add a function to get b=0 TSNR
+                pass
+
         # Get the metadata
         output_file = self.inputs.output_file_name
         image_qc['file_name'] = output_file

qsiprep/workflows/base.py

@@ -22,8 +22,6 @@
 
 from nilearn import __version__ as nilearn_ver
 
-from qsiprep.workflows.fieldmap import pepolar
-
 from ..engine import Workflow
 from ..interfaces import (BIDSDataGrabber, BIDSInfo, BIDSFreeSurferDir,
                           SubjectSummary, AboutSummary, DerivativesDataSink)
@@ -47,7 +45,7 @@ def init_qsiprep_wf(
         subject_list, run_uuid, work_dir, output_dir, bids_dir, ignore, debug,
         low_mem, anat_only, dwi_only, longitudinal, b0_threshold, hires,
         denoise_before_combining, dwi_denoise_window, denoise_method,
-        unringing_method, dwi_no_biascorr, no_b0_harmonization,
+        unringing_method, b1_biascorrect_stage, no_b0_harmonization,
         output_resolution, infant_mode, combine_all_dwis,
         distortion_group_merge, pepolar_method, omp_nthreads, bids_filters,
         force_spatial_normalization, skull_strip_template,
@@ -87,7 +85,7 @@ def init_qsiprep_wf(
                               dwi_denoise_window=7,
                               denoise_method='patch2self',
                               unringing_method='mrdegibbs',
-                              dwi_no_biascorr=False,
+                              b1_biascorrect_stage=False,
                               no_b0_harmonization=False,
                               combine_all_dwis=True,
                               distortion_group_merge='concat',
@@ -154,8 +152,8 @@ def init_qsiprep_wf(
             Either 'dwidenoise', 'patch2self' or 'none'
         unringing_method : str
             algorithm to use for removing Gibbs ringing. Options: none, mrdegibbs
-        dwi_no_biascorr : bool
-            run spatial bias correction (N4) on dwi series
+        b1_biascorrect_stage : str
+            'final', 'none' or 'legacy'
         no_b0_harmonization : bool
             skip rescaling dwi scans to have matching b=0 intensities across scans
         denoise_before_combining : bool
@@ -254,7 +252,7 @@ def init_qsiprep_wf(
             dwi_denoise_window=dwi_denoise_window,
             denoise_method=denoise_method,
             unringing_method=unringing_method,
-            dwi_no_biascorr=dwi_no_biascorr,
+            b1_biascorrect_stage=b1_biascorrect_stage,
             no_b0_harmonization=no_b0_harmonization,
             dwi_only=dwi_only,
             anat_only=anat_only,
@@ -305,7 +303,7 @@ def init_single_subject_wf(
         subject_id, name, reportlets_dir, output_dir, bids_dir, ignore, debug,
         write_local_bvecs, low_mem, dwi_only, anat_only, longitudinal,
         b0_threshold, denoise_before_combining, bids_filters,
-        dwi_denoise_window, denoise_method, unringing_method, dwi_no_biascorr,
+        dwi_denoise_window, denoise_method, unringing_method, b1_biascorrect_stage,
         no_b0_harmonization, infant_mode, combine_all_dwis, raw_image_sdc,
         distortion_group_merge, pepolar_method, omp_nthreads, skull_strip_template,
         force_spatial_normalization, skull_strip_fixed_seed, freesurfer, hires,
@@ -345,7 +343,7 @@ def init_single_subject_wf(
             dwi_denoise_window=7,
             denoise_method='patch2self',
             unringing_method='mrdegibbs',
-            dwi_no_biascorr=False,
+            b1_biascorrect_stage=False,
             no_b0_harmonization=False,
             dwi_only=False,
             anat_only=False,
@@ -407,8 +405,8 @@ def init_single_subject_wf(
             Either 'dwidenoise', 'patch2self' or 'none'
         unringing_method : str
             algorithm to use for removing Gibbs ringing. Options: none, mrdegibbs
-        dwi_no_biascorr : bool
-            run spatial bias correction (N4) on dwi series
+        b1_biascorrect_stage : str
+            'final', 'none' or 'legacy'
         no_b0_harmonization : bool
             skip rescaling dwi scans to have matching b=0 intensities across scans
         denoise_before_combining : bool
@@ -728,7 +726,7 @@ def init_single_subject_wf(
             dwi_denoise_window=dwi_denoise_window,
             denoise_method=denoise_method,
             unringing_method=unringing_method,
-            dwi_no_biascorr=dwi_no_biascorr,
+            b1_biascorrect_stage=b1_biascorrect_stage,
             no_b0_harmonization=no_b0_harmonization,
             denoise_before_combining=denoise_before_combining,
             motion_corr_to=motion_corr_to,
@@ -757,7 +755,6 @@ def init_single_subject_wf(
             name=dwi_preproc_wf.name.replace('dwi_preproc', 'dwi_finalize'),
             output_prefix=output_fname,
             layout=layout,
-            ignore=ignore,
             hmc_model=hmc_model,
             shoreline_iters=shoreline_iters,
             write_local_bvecs=write_local_bvecs,
@@ -767,9 +764,8 @@ def init_single_subject_wf(
             output_resolution=output_resolution,
             output_dir=output_dir,
             omp_nthreads=omp_nthreads,
-            use_syn=use_syn,
-            low_mem=low_mem,
-            pepolar_method=pepolar_method,
+            do_biascorr=b1_biascorrect_stage=='final',
+            b0_threshold=b0_threshold,
             make_intramodal_template=make_intramodal_template,
             source_file=source_file,
             write_derivatives=not merging_distortion_groups)