Re-created bip repository

bip/pipelines/dMRI_diff/dMRI_diff.py

+#!/usr/bin/env python
+#
+# dMRI_diff.py - Pipeline with the dfMRI processing.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import logging
+from bip.utils import redirect_logging
+from bip.pipelines.dMRI_diff import diff_eddy
+from bip.pipelines.dMRI_diff import diff_dtifit
+from bip.pipelines.dMRI_diff import diff_noddi
+from bip.pipelines.dMRI_diff import diff_tbss
+from bip.pipelines.dMRI_diff import diff_bedpostx
+from bip.pipelines.dMRI_diff import diff_autoptx
+
+log = logging.getLogger(__name__)
+
+def add_to_pipeline(ctx, pipe, tree, targets):
+
+    logs_dir=tree.get('logs_dir')
+
+    with redirect_logging('pipe_dMRI_fieldmap', outdir=logs_dir):
+
+        # TODO: Make this configurable
+        cuda_dict_1 = dict(queue="short.qg@@win-uk-biobank-gpu",
+                         coprocessor="cuda",
+                         coprocessor_class="P",
+                         coprocessor_toolkit="8.0")
+        cuda_dict_2 = dict(queue="short.qg@@win-uk-biobank-gpu",
+                         coprocessor="cuda",
+                         coprocessor_class="P",
+                         coprocessor_toolkit="6.5")
+        #pipe(diff_eddy.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        pipe(diff_eddy.run, submit=cuda_dict_1, kwargs={'ctx' : ctx})
+        targets.append('eddy_data')
+        pipe(diff_dtifit.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        targets.append('FA')
+        pipe(diff_noddi.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        targets.append('ISOVF')
+        pipe(diff_tbss.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        targets.append('JHUrois_FA')
+        pipe(diff_bedpostx.run, submit=cuda_dict_2, kwargs={'ctx' : ctx})
+        targets.append('bedpostx_eye')
+        pipe(diff_autoptx.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        targets.append('tractsNorm')
+
+    return pipe, targets

bip/pipelines/dMRI_diff/diff_autoptx.py

+#!/usr/bin/env python
+#
+# diff_autoptx.py - Sub-pipeline with FSL's autoPtx processing of the dMRI.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+# pylint: disable=W0613
+#
+
+import os
+import random
+import logging
+from shutil import copyfile
+from fsl import wrappers
+from bip.utils import redirect_logging
+from pipe_tree import In, Out, Ref, Var
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        autoptx_tract:             Var,
+        bedpostx_eye:              In,
+        TBSS_FA_to_MNI_warp:       In,
+        TBSS_MNI_to_dti_FA_warp:   In,
+        bedpostx_nodif_brain_mask: In,
+        logs_dir:                  Ref,
+        aptx_stop:                 Ref,
+        bedpostx_merged:           Ref,
+        aptx_tracts_sub_dir:       Ref,
+        aptx_tractsinv_sub_dir:    Ref,
+        density_inv:               Ref,
+        waytotal_inv:              Ref,
+        fwdensity:                 Ref,
+        aptx_seed:                 Out,
+        aptx_target:               Out,
+        aptx_exclude:              Out,
+        density:                   Out,
+        waytotal:                  Out,
+        tractsNorm:                Out):
+
+    with redirect_logging('diff_autoptx_' + autoptx_tract.value,
+                          outdir=logs_dir):
+
+        t_name  = autoptx_tract.value
+        #TODO: Allow for customization for this multiplication factor
+        t_seeds = int(ctx.tract_struct[t_name]["num_seeds"] * 300)
+
+        orig_tract_dir = ctx.get_data('dMRI/autoptx/protocols/' + t_name + '/')
+
+        # Does the protocol defines a second
+        # run with inverted seed / target masks?
+        if os.path.exists(orig_tract_dir + 'invert'):
+            symtrack=True
+            if os.path.exists(waytotal):
+                os.remove(waytotal)
+        else:
+            symtrack=False
+
+        # Copy the seed, target and exclude files for the tract
+        copyfile(orig_tract_dir + 'seed.nii.gz',    aptx_seed)
+        copyfile(orig_tract_dir + 'target.nii.gz',  aptx_target)
+        copyfile(orig_tract_dir + 'exclude.nii.gz', aptx_exclude)
+
+        if os.path.exists(waytotal):
+            os.remove(waytotal)
+
+        # Building arguments dictionary for probtrackx
+        kwargs = {
+            'samples' :   bedpostx_merged,
+            'mask' :      bedpostx_nodif_brain_mask,
+            'out' :       'density',
+            'nsamples' :  t_seeds,
+            'avoid' :     aptx_exclude,
+            'opd' :       True,
+            'loopcheck' : True,
+            'forcedir' :  True,
+            'sampvox' :   True,
+            'xfm' :       TBSS_MNI_to_dti_FA_warp,
+            'invxfm' :    TBSS_FA_to_MNI_warp,
+            'rseed' :     int(random.randrange(32767))
+        }
+
+        # Is there a stop criterion defined
+        # in the protocol for this struct?
+        if os.path.exists(orig_tract_dir + 'stop.nii.gz'):
+            kwargs['stop'] = aptx_stop
+            copyfile(orig_tract_dir + 'stop.nii.gz', aptx_stop)
+
+        # Process structure
+        wrappers.probtrackx(seed=aptx_seed,
+                            waypoints=aptx_target,
+                            dir=aptx_tracts_sub_dir,
+                            **kwargs)
+
+        with open(waytotal, 'r', encoding="utf-8") as f:
+            way = int(f.read())
+
+        # Calculate inverse tractography for this tract and
+        # merge runs for forward and inverted tractographies
+        if symtrack:
+            wrappers.probtrackx(seed=aptx_target,
+                                waypoints=aptx_seed,
+                                dir=aptx_tractsinv_sub_dir,
+                                **kwargs)
+
+            copyfile(density, fwdensity)
+            wrappers.fslmaths(density).add(density_inv).run(density)
+
+            with open(waytotal, 'r', encoding="utf-8") as f:
+                way1 = int(f.read())
+            with open(waytotal_inv, 'r', encoding="utf-8") as f:
+                way2 = int(f.read())
+            if os.path.exists(waytotal):
+                os.remove(waytotal)
+            with open(waytotal, 'wt', encoding="utf-8") as f:
+                way = way1 + way2
+                f.write(f'{way}\n')
+
+        wrappers.fslmaths(density).div(way).range().run(tractsNorm)

bip/pipelines/dMRI_diff/diff_bedpostx.py

+#!/usr/bin/env python
+#
+# diff_bedpostx.py - Sub-pipeline with FSL's bedpostx processing of the dMRI.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+# pylint: disable=W0613
+#
+
+import os
+import logging
+from shutil import copyfile
+from fsl import wrappers
+from bip.utils import redirect_logging
+from pipe_tree import In, Out, Ref
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        AP_bval:                   In,
+        AP_bvec:                   In,
+        eddy_data_ud:              In,
+        eddy_nodif:                In,
+        eddy_nodif_brain_mask:     In,
+        logs_dir:                  Ref,
+        eddy_dir:                  Ref,
+        bedpostx_dir:              Ref,
+        bedpostx_log_part:         Ref,
+        bedpostx_data_0:           Ref,
+        bedpostx_bvals:            Out,
+        bedpostx_bvecs:            Out,
+        bedpostx_nodif_brain:      Out,
+        bedpostx_nodif_brain_mask: Out,
+        bedpostx_nvox:             Out,
+        bedpostx_eye:              Out,
+        bp_logs_gpu_dir:           Out):
+
+    with redirect_logging('diff_bedpostx', outdir=logs_dir):
+
+        w = os.getcwd() + "/"
+
+        # Copy of bvals and bvecs
+        copyfile(src=AP_bval, dst=bedpostx_bvals)
+        copyfile(src=AP_bvec, dst=bedpostx_bvecs)
+        copyfile(src=eddy_nodif_brain_mask, dst=bedpostx_nodif_brain_mask)
+
+        wrappers.fslmaths(eddy_nodif).mas(eddy_nodif_brain_mask).run(bedpostx_nodif_brain)
+
+        nvox=wrappers.fslstats(bedpostx_nodif_brain_mask).V.run()[0]
+        with open(bedpostx_nvox, 'wt', encoding="utf-8") as f:
+            f.write(str(nvox))
+
+        wrappers.split_parts_gpu(Datafile=w + eddy_data_ud,
+                                 Maskfile=w + eddy_nodif_brain_mask,
+                                 Bvalsfile=w + bedpostx_bvals,
+                                 Bvecsfile=w + bedpostx_bvecs,
+                                 TotalNumParts="1",
+                                 OutputDirectory=w + bedpostx_dir)
+
+        wrappers.xfibres_gpu(data=w + bedpostx_data_0,
+                             mask=w + bedpostx_nodif_brain_mask,
+                             bvals=w + bedpostx_bvals, bvecs=w + bedpostx_bvecs,
+                             forcedir=True, logdir=w + bedpostx_log_part,
+                             nfibres="3", fudge="1", burnin="3000",
+                             njumps="1250", sampleevery="25", model="2",
+                             cnonlinear=True, subjdir=w + eddy_dir,
+                             idpart="0", nparts="1",
+                             numvoxels=str(nvox))
+
+        wrappers.bedpostx_postproc_gpu(data=w + bedpostx_data_0,
+                                       mask=w + bedpostx_nodif_brain_mask,
+                                       bvals=w + bedpostx_bvals,
+                                       bvecs=w + bedpostx_bvecs, forcedir=True,
+                                       logdir=w + bedpostx_log_part, nf="3",
+                                       fudge="1", bi="3000", nj="1250", se="25",
+                                       model="2", cnonlinear=True,
+                                       TotalNumVoxels=str(nvox),
+                                       SubjectDir=w + eddy_dir,
+                                       TotalNumParts="1", bindir=ctx.FSLDIR)
+
+        if os.path.exists(bedpostx_data_0):
+            os.remove(bedpostx_data_0)

bip/pipelines/dMRI_diff/diff_dtifit.py

+#!/usr/bin/env python
+#
+# diff_dtifit.py - Sub-pipeline with FSL's dtifit processing of the dMRI.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+# pylint: disable=W0613
+#
+
+import logging
+from shutil import copyfile
+import numpy as np
+from fsl import wrappers
+from bip.utils import redirect_logging
+from pipe_tree import In, Out, Ref
+from gradunwarp.core.gradient_unwarp_apply import gradient_unwarp_apply
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        eddy_data:                 In,
+        eddy_bvals:                In,
+        eddy_bvecs:                In,
+        eddy_nodif_brain_mask_ud:  In,
+        logs_dir:                  Ref,
+        eddy_data_GDC:             Ref,
+        dtifit_output_prefix:      Ref,
+        eddy_data_ud:              Out,
+        eddy_data_ud_warp:         Out,
+        eddy_data_ud_1_shell:      Out,
+        eddy_data_ud_1_shell_bval: Out,
+        eddy_data_ud_1_shell_bvec: Out,
+        FA:                        Out):
+
+    with redirect_logging('diff_dtifit', outdir=logs_dir):
+
+
+        if ctx.gdc != '':
+            #Calculate and apply the Gradient Distortion Unwarp
+            # TODO: Review the "half=True" in next version
+            gradient_unwarp_apply(WD=eddy_data_GDC,
+                                  infile=eddy_data,
+                                  outfile=eddy_data_ud,
+                                  owarp=eddy_data_ud_warp,
+                                  gradcoeff=ctx.gdc,
+                                  vendor='siemens', nojac=True, half=True)
+        else:
+            copyfile(src=eddy_data, dst=eddy_data_ud)
+
+
+        #Correct input for dtifit using one shell
+        list_vals = []
+        b_vals = np.loadtxt(eddy_bvals)
+        b_vecs = np.loadtxt(eddy_bvecs)
+        lim_file = ctx.get_data('dMRI/b0_threshold.txt')
+
+        lim = int(np.loadtxt(lim_file))
+
+        # TODO: This can be modifiable by the user
+        b_value_shell_to_keep = 1000
+
+        for count, b_val in enumerate(b_vals):
+            # These are the b0 images: They are automatically included
+            if b_val < lim:
+                list_vals.append(count)
+            else:
+                if abs(b_val - b_value_shell_to_keep) < lim:
+                    list_vals.append(count)
+
+        new_b_vals = b_vals[list_vals]
+        new_b_vecs = b_vecs[:,list_vals]
+
+        np.savetxt(eddy_data_ud_1_shell_bval, new_b_vals, newline=" ", fmt='%.1f')
+        np.savetxt(eddy_data_ud_1_shell_bvec, new_b_vecs, fmt='%.6f')
+
+        wrappers.fslselectvols(src=eddy_data_ud, out=eddy_data_ud_1_shell,
+                               vols=list_vals)
+
+        wrappers.dtifit(data  = eddy_data_ud_1_shell,
+                        bvecs = eddy_data_ud_1_shell_bvec,
+                        bvals = eddy_data_ud_1_shell_bval,
+                        mask  = eddy_nodif_brain_mask_ud,
+                        out   = dtifit_output_prefix,
+                        save_tensor = True)

bip/pipelines/dMRI_diff/diff_eddy.py

+#!/usr/bin/env python
+#
+# diff_eddy.py - Sub-pipeline with FSL's eddy processing of the dMRI.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+# pylint: disable=W0613
+#
+
+import os
+import logging
+from shutil import copyfile
+import nibabel as nib
+from fsl import wrappers
+from bip.utils import redirect_logging
+from pipe_tree import In, Out, Ref
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        AP:                             In,
+        AP_bval:                        In,
+        AP_bvec:                        In,
+        AP_best_index:                  In,
+        fieldmap_iout_mean:             In,
+        fieldmap_mask:                  In,
+        fieldmap_mask_ud:               In,
+        acqparams:                      In,
+        logs_dir:                       Ref,
+        eddy_data_prefix:               Ref,
+        fieldmap_out_prefix:            Ref,
+        eddy_AP:                        Out,
+        eddy_nodif:                     Out,
+        eddy_nodif_brain_mask:          Out,
+        eddy_nodif_brain_mask_ud:       Out,
+        eddy_bvals:                     Out,
+        eddy_bvecs:                     Out,
+        eddy_index:                     Out,
+        eddy_data:                      Out):
+
+    with redirect_logging('diff_eddy', outdir=logs_dir):
+
+        # Creates links
+        # TODO: These links are NOT relative. This may cause future problems.
+        if not os.path.exists(eddy_AP):
+            os.symlink(src="../../../" + AP, dst=eddy_AP)
+        if not os.path.exists(eddy_nodif):
+            os.symlink(src="../../../" + fieldmap_iout_mean, dst=eddy_nodif)
+        if not os.path.exists(eddy_nodif_brain_mask):
+            os.symlink(src="../../../" + fieldmap_mask, dst=eddy_nodif_brain_mask)
+        if not os.path.exists(eddy_nodif_brain_mask_ud):
+            os.symlink(src="../../../" + fieldmap_mask_ud, dst=eddy_nodif_brain_mask_ud)
+
+        # Generation of FSF file
+        copyfile(src=AP_bval, dst=eddy_bvals)
+        copyfile(src=AP_bvec, dst=eddy_bvecs)
+
+        AP_img    = nib.load(eddy_AP)
+        AP_numvol = AP_img.header['dim'][4]
+
+        indices = ("1 " * AP_numvol).strip() + "\n"
+
+        with open(AP_best_index, 'r', encoding="utf-8") as f:
+            AP_index = int(f.read())
+
+        with open(eddy_index, 'wt', encoding="utf-8") as f:
+            f.write(indices)
+
+        wrappers.eddy(imain=eddy_AP,
+                      mask=eddy_nodif_brain_mask,
+                      topup=fieldmap_out_prefix,
+                      acqp=acqparams,
+                      index=eddy_index,
+                      bvecs=eddy_bvecs,
+                      bvals=eddy_bvals,
+                      out=eddy_data_prefix,
+                      ref_scan_no=AP_index,
+                      flm="quadratic",
+                      resamp="jac",
+                      slm="linear",
+                      niter=8,
+                      fwhm=[10,8,4,2,0,0,0,0],
+                      ff=10,
+                      sep_offs_move=True,
+                      nvoxhp=1000,
+                      very_verbose=True,
+                      repol=True,
+                      rms=True)

bip/pipelines/dMRI_diff/diff_noddi.py

+#!/usr/bin/env python
+#
+# diff_noddi.py - Sub-pipeline with AMICO processing of the dMRI.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import os
+import logging
+from shutil import move, rmtree
+import zipfile
+import gzip
+import shutil
+import amico
+from bip.utils import redirect_logging, tempdir
+from pipe_tree import In, Out, Ref
+
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        eddy_data_ud:              In,
+        eddy_nodif_brain_mask_ud:  In,
+        eddy_bvals:                In,
+        eddy_bvecs:                In,
+        logs_dir:                  Ref,
+        tmp_dir:                   Ref,
+        NODDI_dir:                 Ref,
+        AMICO_dir:                 Ref,
+        NODDI_kernels_dir:         Ref,
+        NODDI_scheme:              Out,
+        ICVF:                      Out,
+        ISOVF:                     Out,
+        OD:                        Out,
+        NODDI_dir_file:            Out):
+
+    with redirect_logging('diff_noddi', outdir=logs_dir),\
+         tempdir(tmp_dir):
+        amtmp      = tmp_dir + '/amtmp.nii'
+        amtmp_mask = tmp_dir + '/amtmp_mask.nii'
+
+        # In case the kernels were already calculated for this study:
+        # TODO: Verify that this is regular UKB data.
+        # Otherwise, this should not be done
+        if not os.path.exists('./kernels/'):
+            if os.path.exists(ctx.get_data('dMRI/kernels.zip')):
+                with zipfile.ZipFile(ctx.get_data('dMRI/kernels.zip'), 'r') as zip_ref:
+                    zip_ref.extractall("./")
+
+        with gzip.open(eddy_data_ud, 'rb') as f_in:
+            with open(amtmp, 'wb') as f_out:
+                shutil.copyfileobj(f_in, f_out)
+
+        with gzip.open(eddy_nodif_brain_mask_ud, 'rb') as f_in:
+            with open(amtmp_mask , 'wb') as f_out:
+                shutil.copyfileobj(f_in, f_out)
+
+        amico.setup()
+        amico.util.fsl2scheme(bvalsFilename=eddy_bvals,
+                              bvecsFilename=eddy_bvecs,
+                              schemeFilename=NODDI_scheme,
+                              bStep=50)
+
+        ae = amico.Evaluation("./", "./", output_path=NODDI_dir)
+        ae.load_data(amtmp, NODDI_scheme, mask_filename=amtmp_mask,
+                     b0_min_signal=1e-4)
+        ae.set_model("NODDI")
+        ae.generate_kernels(ndirs=500,regenerate=True)
+        ae.load_kernels()
+        ae.fit()
+        ae.save_results()
+
+        files = [ICVF, ISOVF, OD, NODDI_dir_file]
+
+        # Move files to proper place
+        for file_n in files:
+            old_name = file_n.replace("/NODDI_","/FIT_")
+            move(old_name, file_n)
+
+        # Delete created folders
+        for dir_to_delete in [AMICO_dir, NODDI_kernels_dir]:
+            if os.path.exists(dir_to_delete):
+                rmtree(dir_to_delete)

bip/pipelines/dMRI_diff/diff_tbss.py

+#!/usr/bin/env python
+#
+# diff_tbss.py - Sub-pipeline with FSL's TBSS processing of the dMRI.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+# pylint: disable=W0613,C0206
+#
+
+import os
+import logging
+from shutil import copyfile
+import nibabel as nib
+from fsl import wrappers
+from bip.utils import redirect_logging, tempdir
+from pipe_tree import In, Out, Ref
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        FA:                          In,
+        ISOVF:                       In(optional=True),
+        logs_dir:                    Ref,
+        tmp_dir:                     Ref,
+        TBSS_MNI_to_dti_FA_warp_msf: Ref,
+        dtifit_prefix:               Ref,
+        NODDI_prefix:                Ref,
+        TBSS_prefix:                 Ref,
+        JHUrois_prefix:              Ref,
+        TBSS_FA_to_MNI_int:          Ref,
+        TBSS_FA_dir_FA:              Out,
+        TBSS_FA_dir_FA_mask:         Out,
+        TBSS_FA:                     Out,
+        TBSS_MNI:                    Out,
+        TBSS_FA_to_MNI_affine_mat:   Out,
+        TBSS_FA_to_MNI_warp:         Out,
+        TBSS_FA_to_MNI_int_txt:      Out,
+        TBSS_log_1:                  Out,
+        TBSS_log_2:                  Out,
+        TBSS_log_3:                  Out,
+        TBSS_dti_FA_to_MNI:          Out,
+        TBSS_MNI_to_dti_FA_warp:     Out,
+        TBSS_all_FA:                 Out,
+        TBSS_mean_FA_mask:           Out,
+        TBSS_mean_FA:                Out,
+        TBSS_mean_FA_skeleton:       Out,
+        TBSS_mean_FA_skeleton_mask:  Out,
+        TBSS_all_FA_skeletonised:    Out,
+        JHUrois_FA:                  Out):
+
+    with redirect_logging('diff_tbss', outdir=logs_dir),\
+         tempdir(tmp_dir):
+
+        TBSS_FA_tmp = tmp_dir + '/TBSS_FA_tmp.nii.gz'
+        TBSS_FA_to_MNI_warp_s1 = tmp_dir + '/FA_to_MNI_warp_s1.nii.gz'
+        TBSS_FA_to_MNI_warp_s2 = tmp_dir + '/FA_to_MNI_warp_s2.nii.gz'
+
+        # Creates links
+        # TODO: These links are NOT relative. This may cause future problems.
+        if not os.path.exists(TBSS_FA):
+            os.symlink(src="../../../" + FA, dst=TBSS_FA)
+        if not os.path.exists(TBSS_MNI):
+            os.symlink(src=ctx.get_standard('FMRIB58_FA_1mm.nii.gz'), dst=TBSS_MNI)
+        FA_img     = nib.load(TBSS_FA)
+        x, y, z = FA_img.header['dim'][1:4]
+
+        # erode a little and zero end slices
+        wrappers.fslmaths(FA).min(1).ero().roi(1,x,1,y,1,z,0,1).run(TBSS_FA_dir_FA)
+
+        # create mask (for use in FLIRT & FNIRT)
+        wrappers.fslmaths(TBSS_FA_dir_FA).bin().run(TBSS_FA_dir_FA_mask)
+        #TODO: investigate how to add the -odt=char option
+        wrappers.fslmaths(TBSS_FA_dir_FA_mask).dilD(2).sub(1).abs().add(TBSS_FA_dir_FA_mask).run(TBSS_FA_dir_FA_mask)
+
+        if not os.path.exists(TBSS_MNI_to_dti_FA_warp_msf):
+            #New Optimal Registration
+            wrappers.flirt(ref=TBSS_MNI,
+                           src=TBSS_FA_dir_FA,
+                           inweight=TBSS_FA_dir_FA_mask,
+                           omat=TBSS_FA_to_MNI_affine_mat)
+
+            # perform FNIRT cascade of registrations
+            wc1 = ctx.get_data('dMRI/dMRI_reg_optimal_parameters/oxford_s1.cnf')
+            wc2 = ctx.get_data('dMRI/dMRI_reg_optimal_parameters/oxford_s2.cnf')
+            wc3 = ctx.get_data('dMRI/dMRI_reg_optimal_parameters/oxford_s3.cnf')
+
+            wrappers.fnirt(ref=TBSS_MNI, src=TBSS_FA_dir_FA,
+                           cout=TBSS_FA_to_MNI_warp_s1, config=wc1,
+                           aff=TBSS_FA_to_MNI_affine_mat,
+                           intout=TBSS_FA_to_MNI_int, logout=TBSS_log_1)
+            wrappers.fnirt(ref=TBSS_MNI, src=TBSS_FA_dir_FA,
+                           cout=TBSS_FA_to_MNI_warp_s2, config=wc2,
+                           inwarp=TBSS_FA_to_MNI_warp_s1,
+                           intin=TBSS_FA_to_MNI_int_txt, logout=TBSS_log_2)
+            wrappers.fnirt(ref=TBSS_MNI, src=TBSS_FA_dir_FA,
+                           cout=TBSS_FA_to_MNI_warp, config=wc3,
+                           iout=TBSS_dti_FA_to_MNI,
+                           inwarp=TBSS_FA_to_MNI_warp_s2,
+                           intin=TBSS_FA_to_MNI_int_txt, logout=TBSS_log_3)
+            wrappers.invwarp(ref=TBSS_FA_dir_FA, warp=TBSS_FA_to_MNI_warp,
+                             out=TBSS_MNI_to_dti_FA_warp)
+
+
+            # now estimate the mean deformation
+            wrappers.fslmaths(TBSS_FA_to_MNI_warp).sqr().Tmean().run(TBSS_FA_tmp)
+            mean    = wrappers.fslstats(TBSS_FA_tmp).M.run()
+            median  = wrappers.fslstats(TBSS_FA_tmp).P(50).run()
+
+            with open(TBSS_MNI_to_dti_FA_warp_msf, 'wt', encoding="utf-8") as f:
+                f.write(f'{mean} {median}')
+
+        wrappers.applywarp(src=TBSS_FA_dir_FA, out=TBSS_dti_FA_to_MNI,
+                           ref=TBSS_MNI, warp=TBSS_FA_to_MNI_warp,
+                           rel=True, interp='trilinear')
+
+        copyfile(src=TBSS_dti_FA_to_MNI, dst=TBSS_all_FA)
+
+        # Create mean FA
+        skel = ctx.get_standard('FMRIB58_FA-skeleton_1mm.nii.gz')
+        wrappers.fslmaths(TBSS_all_FA).bin().mul(TBSS_MNI).bin().run(TBSS_mean_FA_mask)
+        wrappers.fslmaths(TBSS_all_FA).mas(TBSS_mean_FA_mask).run(TBSS_all_FA)
+        wrappers.fslmaths(TBSS_MNI).mas(TBSS_mean_FA_mask).run(TBSS_mean_FA)
+        wrappers.fslmaths(skel).mas(TBSS_mean_FA_mask).run(TBSS_mean_FA_skeleton)
+
+        # Create Skeleton mask
+        # TODO: This threshold could be configurable by user
+        thresh = 2000
+        wrappers.fslmaths(TBSS_mean_FA_skeleton).thr(thresh).bin().run(TBSS_mean_FA_skeleton_mask)
+        wrappers.fslmaths(TBSS_all_FA).mas(TBSS_mean_FA_skeleton_mask).run(TBSS_all_FA_skeletonised)
+
+        atlas = ctx.get_atlas('JHU/JHU-ICBM-labels-1mm.nii.gz')
+        mean = wrappers.fslstats(TBSS_all_FA_skeletonised, K=atlas).M.run()
+
+        with open(JHUrois_FA, 'wt', encoding="utf-8") as f:
+            f.write(f'{mean}')
+
+        # Applying to the outputs of both dtifit and noddi
+        dict_suffix={dtifit_prefix: ['L1', 'L2', 'L3', 'MO', 'MD'],
+                     NODDI_prefix: ['ICVF', 'OD', 'ISOVF']}
+
+        for prefix in dict_suffix:
+
+            for suffix in dict_suffix[prefix]:
+
+                d_input = prefix + suffix + ".nii.gz"
+
+                if os.path.exists(d_input):
+                    d_output      = TBSS_prefix + suffix +".nii.gz"
+                    d_output_skel = TBSS_prefix + suffix +"_skeletonised.nii.gz"
+                    d_output_txt  = JHUrois_prefix + suffix + '.txt'
+                    wrappers.applywarp(src=d_input, out=d_output, ref=TBSS_MNI,
+                                       warp=TBSS_FA_to_MNI_warp, rel=True,
+                                       interp='trilinear')
+                    wrappers.fslmaths(d_output).mas(TBSS_mean_FA_skeleton_mask).run(d_output_skel)
+                    mean = wrappers.fslstats( d_output_skel, K=atlas).M.run()
+
+                    with open(d_output_txt, 'wt', encoding="utf-8") as f:
+                        f.write(f'{mean}')

bip/pipelines/dMRI_fieldmap/dMRI_fieldmap.py

+#!/usr/bin/env python
+#
+# dMRI_fieldmap.py - Pipeline with the fieldmap processing.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import logging
+from bip.utils import redirect_logging
+from bip.pipelines.dMRI_fieldmap import fieldmap_pre_topup
+from bip.pipelines.dMRI_fieldmap import fieldmap_topup
+from bip.pipelines.dMRI_fieldmap import fieldmap_post_topup
+
+log = logging.getLogger(__name__)
+
+def add_to_pipeline(ctx, pipe, tree, targets):
+
+    logs_dir=tree.get('logs_dir')
+
+    with redirect_logging('pipe_dMRI_fieldmap', outdir=logs_dir):
+        pipe(fieldmap_pre_topup.run, submit=dict(jobtime=200),
+             kwargs={'ctx' : ctx})
+        targets.append('total_B0_PA')
+        pipe(fieldmap_topup.run, submit=dict(jobtime=200),
+             kwargs={'ctx' : ctx})
+        targets.append('fieldmap_fout')
+        pipe(fieldmap_post_topup.run, submit=dict(jobtime=200),
+             kwargs={'ctx' : ctx})
+        targets.append('fieldmap_mask')
+
+    return pipe, targets

bip/pipelines/dMRI_fieldmap/fieldmap_post_topup.py

+#!/usr/bin/env python
+#
+# fieldmap_post_topup.py - Sub-pipeline with the fieldmap processing after topup
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import logging
+from shutil import copyfile
+import numpy as np
+import nibabel as nib
+from fsl import wrappers
+from bip.utils import redirect_logging, tempdir
+from pipe_tree import In, Out, Ref
+from gradunwarp.core.gradient_unwarp_apply import gradient_unwarp_apply
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        T1:                             In,
+        T1_brain:                       In,
+        T1_brain_mask:                  In,
+        T1_fast_WM_mask:                In,
+        T1_brain_mask_dil:              Out,
+        B0_AP:                          In,
+        B0_PA:                          In,
+        acqparams:                      In,
+        logs_dir:                       Ref,
+        tmp_dir:                        Ref,
+        fieldmap_out_prefix:            Ref,
+        fieldmap_GDC:                   Ref,
+        fieldmap_fout:                  Ref,
+        T1_to_fieldmap_iout_mat:        Out,
+        fieldmap_iout:                  Out,
+        fieldmap_iout_mean:             Out,
+        fieldmap_iout_mean_ud:          Out,
+        fieldmap_iout_mean_ud_warp:     Out,
+        fieldmap_iout_to_T1:            Out,
+        fieldmap_fout_to_T1:            Out,
+        fieldmap_fout_to_T1_brain:      Out,
+        fieldmap_fout_to_T1_brain_rad:  Out,
+        fieldmap_iout_to_T1_mat:        Out,
+        fieldmap_mask:                  Out,
+        fieldmap_mask_ud:               Out,
+        fieldmap_iout_mean_ud_inv_warp: Out):
+
+    with redirect_logging('fieldmap_post_topup', outdir=logs_dir),\
+               tempdir(tmp_dir):
+
+        B0_AP_corr_tmp  = tmp_dir + '/B0_AP_corr_tmp.nii.gz'
+        B0_PA_corr_tmp  = tmp_dir + '/B0_PA_corr_tmp.nii.gz'
+        B0_AP_fixed_tmp = tmp_dir + '/B0_AP_fixed_tmp.nii.gz'
+        B0_PA_fixed_tmp = tmp_dir + '/B0_PA_fixed_tmp.nii.gz'
+        fieldmap_tmp    = tmp_dir + '/fieldmap_tmp.nii.gz'
+
+        #Get corrected B=0 for AP and PA and average them avoiding zero values
+        wrappers.applytopup(imain=B0_AP, datain=acqparams, index=1,
+                            topup=fieldmap_out_prefix,
+                            out=B0_AP_corr_tmp, method="jac")
+        wrappers.applytopup(imain=B0_PA, datain=acqparams, index=2,
+                            topup=fieldmap_out_prefix,
+                            out=B0_PA_corr_tmp, method="jac")
+
+        B0_AP_corr_tmp_img = nib.load(B0_AP_corr_tmp)
+        B0_PA_corr_tmp_img = nib.load(B0_PA_corr_tmp)
+
+        B0_AP_corr_tmp_imgf = B0_AP_corr_tmp_img.get_fdata()
+        B0_PA_corr_tmp_imgf = B0_PA_corr_tmp_img.get_fdata()
+
+        #Get a list of indices the zero-valued voxels for AP and PA
+        x1,y1,z1 = np.where(B0_AP_corr_tmp_imgf==0)
+        x2,y2,z2 = np.where(B0_PA_corr_tmp_imgf==0)
+
+        #For AP voxels with 0 value,, get the values in PA  (And viceversa)
+        B0_AP_fixed_tmp_imgf = B0_AP_corr_tmp_imgf
+        B0_PA_fixed_tmp_imgf = B0_PA_corr_tmp_imgf
+        B0_AP_fixed_tmp_imgf[x1,y1,z1] = B0_PA_corr_tmp_imgf[x1,y1,z1]
+        B0_PA_fixed_tmp_imgf[x2,y2,z2] = B0_AP_corr_tmp_imgf[x2,y2,z2]
+
+        B0_AP_fixed_tmp_img = nib.Nifti1Image(B0_AP_fixed_tmp_imgf,
+                                              B0_AP_corr_tmp_img.affine,
+                                              B0_AP_corr_tmp_img.header)
+
+        B0_PA_fixed_tmp_img = nib.Nifti1Image(B0_PA_fixed_tmp_imgf,
+                                              B0_PA_corr_tmp_img.affine,
+                                              B0_PA_corr_tmp_img.header)
+
+        nib.save(B0_AP_fixed_tmp_img, B0_AP_fixed_tmp)
+        nib.save(B0_PA_fixed_tmp_img, B0_PA_fixed_tmp)
+
+        #Merge and average them
+        wrappers.fslmerge("t", fieldmap_iout, B0_AP_fixed_tmp, B0_PA_fixed_tmp)
+        wrappers.fslmaths(fieldmap_iout).Tmean().run(fieldmap_iout_mean)
+
+
+        if ctx.gdc != '':
+            #Calculate and apply the Gradient Distortion Unwarp
+            # TODO: Review the "half=True" in next version
+            gradient_unwarp_apply(WD=fieldmap_GDC, infile=fieldmap_iout_mean,
+                                  outfile=fieldmap_iout_mean_ud,
+                                  owarp=fieldmap_iout_mean_ud_warp,
+                                  gradcoeff=ctx.gdc,
+                                  vendor='siemens', nojac=True, half=True)
+        else:
+            copyfile(src=fieldmap_iout_mean, dst=fieldmap_iout_mean_ud)
+
+        wrappers.fslcpgeom(B0_AP,fieldmap_fout)
+
+        # Get the topup iout (magnitude) to struct space and
+        # apply the transformation to fout (fieldmap)
+        wrappers.epi_reg(epi=fieldmap_iout_mean_ud, t1=T1, t1brain=T1_brain,
+                         out=fieldmap_iout_to_T1, wmseg=T1_fast_WM_mask)
+        wrappers.applywarp(src=fieldmap_fout, ref=T1,
+                           out=fieldmap_fout_to_T1,
+                           w=fieldmap_iout_mean_ud_warp,
+                           postmat=fieldmap_iout_to_T1_mat,
+                           rel=True, interp='spline')
+
+        # Mask the warped fout (fieldmap) using T1 brain mask
+        wrappers.fslmaths(fieldmap_fout_to_T1).mul(T1_brain_mask).run(fieldmap_fout_to_T1_brain)
+
+        # Multiply the warped & masked fout (fieldmap) by 2*Pi to have it in radians
+        wrappers.fslmaths(fieldmap_fout_to_T1_brain).mul(2*np.pi).run(fieldmap_fout_to_T1_brain_rad)
+
+        # Generate a mask for topup output by inverting the
+        # previous registration and applying it to T1 brain mask
+        wrappers.invxfm(inmat=fieldmap_iout_to_T1_mat, omat=T1_to_fieldmap_iout_mat)
+        wrappers.fslmaths(T1_brain_mask).thr(0.1).run(fieldmap_tmp)
+        wrappers.fslmaths(fieldmap_tmp).kernel("sphere",1.1).dilF().run(fieldmap_tmp)
+        wrappers.fslmaths(fieldmap_tmp).bin().fillh().run(T1_brain_mask_dil)
+        wrappers.applyxfm(src=T1_brain_mask_dil, ref=fieldmap_iout_mean,
+                          mat=T1_to_fieldmap_iout_mat, out=fieldmap_mask_ud,
+                          interp="trilinear")
+        wrappers.fslmaths(fieldmap_mask_ud).thr(0.25).bin().run(fieldmap_mask_ud)
+
+        #Warp the dilated T1 brain mask to the Gradient Distorted space by inverting
+        wrappers.invwarp(ref=fieldmap_iout_mean, warp=fieldmap_iout_mean_ud_warp,
+                         out=fieldmap_iout_mean_ud_inv_warp)
+        wrappers.applywarp(src=T1_brain_mask_dil, ref=fieldmap_iout_mean,
+                           out=fieldmap_mask,
+                           premat=T1_to_fieldmap_iout_mat,
+                           w=fieldmap_iout_mean_ud_inv_warp,
+                           rel=True, interp='trilinear')
+        wrappers.fslmaths(fieldmap_mask).thr(0.25).bin().fillh().run(fieldmap_mask)

bip/pipelines/dMRI_fieldmap/fieldmap_pre_topup.py

+#!/usr/bin/env python
+#
+# fieldmap_pre_topup.py - Sub-pipeline with the fieldmap processing before topup
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+# pylint: disable=W0613,R0912
+#
+
+import glob
+import logging
+from shutil import copyfile
+import numpy as np
+import nibabel as nib
+from fsl import wrappers
+from bip.utils import redirect_logging, tempdir
+from bip.commands import bb_get_b0s, get_dwell_time
+from bip.commands.bb_read_json_field import bb_read_json_field
+from pipe_tree import In, Out, Ref
+
+
+log = logging.getLogger(__name__)
+
+def choose_best_B0(ctx, img, bval, total_B0, indices, tmp, best_index, B0,
+                   fieldmap_flag, tmp_dir):
+
+    blah     = tmp_dir + '/blah.nii.gz'
+    blah_mat = tmp_dir + '/blah.mat'
+
+    b0_threshold=int(np.loadtxt(ctx.get_data('dMRI/b0_threshold.txt')))
+
+    bvals = np.loadtxt(bval)
+
+    if type(bvals) is not np.ndarray:
+        num_vols = 1
+    else:
+        num_vols = len(np.where(bvals <= b0_threshold)[0])
+
+    bb_get_b0s.bb_get_b0s(img, bval, total_B0, indices, num_vols,  b0_threshold)
+    ind_values = [int(x) for x in np.loadtxt(indices)]
+
+    wrappers.fslsplit(total_B0, tmp)
+
+    files = glob.glob(tmp + '*')
+    files.sort()
+
+    N = len(files)
+
+    if N == 1:
+        with open(best_index, 'wt', encoding="utf-8") as f:
+            f.write('0\n')
+        copyfile(src=total_B0, dst=B0)
+    else:
+        scores= np.zeros([N,N])
+        for i in range(N):
+            for j in range(N):
+                if j > i:
+                    wrappers.flirt(src=files[i], ref=files[j], nosearch=True,
+                                   dof=6, out=blah, omat=blah_mat)
+                    im1  = nib.load(blah)
+                    im2  = nib.load(files[j])
+                    im1f = im1.get_fdata().flatten()
+                    im2f = im2.get_fdata().flatten()
+                    corr = np.corrcoef(im1f,im2f)[0,1]
+                    scores[i,j] = corr
+                    scores[j,i] = corr
+
+        final_scores = np.zeros(N)
+        for i in range(N):
+            final_scores[i]=np.sum(scores[:,i]) / (N - 1)
+
+        best_ind = final_scores.argmax()
+        best_score = final_scores[best_ind]
+
+        if final_scores[0] > 0.98:
+            best_ind = 0
+            best_score = final_scores[best_ind]
+        if best_score < 0.95:
+            with open(fieldmap_flag, 'wt', encoding="utf-8") as f:
+                f.write(f'The best score was {best_score}.' +
+                 'The subject will probably have problems with this DWI data\n')
+        with open(best_index, 'wt', encoding="utf-8") as f:
+            f.write(f'{ind_values[best_ind]}\n')
+            copyfile(src=files[best_ind], dst=B0)
+
+
+def generate_acqparams(AP, PA, AP_json, PA_json, numAP, numPA, acqparams):
+
+    enc_dir_AP = bb_read_json_field(fileName=AP_json,
+                                 fieldName="PhaseEncodingDirection")
+    enc_dir_PA = bb_read_json_field(fileName=PA_json,
+                                 fieldName="PhaseEncodingDirection")
+
+    # TODO: Write a way to evaluate that they are indeed opposite directions
+
+    # Generating the specific strings for the acqparams file for AP
+    # according to the Phase Encoding Direction in the AP json file
+    if enc_dir_AP == "j-":
+        dim_AP = 2
+        values_AP = "0 -1 0"
+    elif enc_dir_AP == "j":
+        dim_AP = 2
+        values_AP = "0 1 0"
+    elif enc_dir_AP == "i-":
+        dim_AP = 1
+        values_AP = "-1 0 0"
+    elif enc_dir_AP == "i":
+        dim_AP = 1
+        values_AP = "1 0 0"
+    elif enc_dir_AP == "k-":
+        dim_AP = 3
+        values_AP = "0 0 -1"
+    else:
+        dim_AP = 3
+        values_AP = "0 0 1"
+
+    # Generating the specific strings for the acqparams file for PA
+    # according to the Phase Encoding Direction in the PA json file
+    if enc_dir_PA == "j-":
+        #dim_PA = 2
+        values_PA = "0 -1 0"
+    elif enc_dir_PA == "j":
+        #dim_PA = 2
+        values_PA = "0 1 0"
+    elif enc_dir_PA == "i-":
+        #dim_PA = 1
+        values_PA = "-1 0 0"
+    elif enc_dir_PA == "i":
+        #dim_PA = 1
+        values_PA = "1 0 0"
+    elif enc_dir_PA == "k-":
+        #dim_PA = 3
+        values_PA = "0 0 -1"
+    else:
+        #dim_PA = 3
+        values_PA = "0 0 1"
+
+
+    imAP = nib.load(AP)
+    numlines = imAP.header['dim'][dim_AP]
+
+    dtiDwell = get_dwell_time.get_dt(AP, AP_json)
+
+    topupValue = (dtiDwell * (numlines -1)) / 1000.0
+
+    values_AP = values_AP + " " + str(topupValue) + "\n"
+    values_PA = values_PA + " " + str(topupValue) + "\n"
+
+    with open(acqparams, 'wt', encoding="utf-8") as f:
+        for i in range(numAP):
+            f.write(values_AP)
+        for i in range(numPA):
+            f.write(values_PA)
+
+
+def run(ctx,
+        AP:               In,
+        PA:               In,
+        AP_bval:          In,
+        PA_bval:          In,
+        AP_json:          In,
+        PA_json:          In,
+        AP_indices:       Out,
+        PA_indices:       Out,
+        total_B0_AP:      Out,
+        total_B0_PA:      Out,
+        B0_AP:            Out,
+        B0_PA:            Out,
+        AP_best_index:    Out,
+        PA_best_index:    Out,
+        B0_AP_PA:         Out,
+        acqparams:        Out,
+        AP_fieldmap_flag: Ref,
+        PA_fieldmap_flag: Ref,
+        logs_dir:         Ref,
+        tmp_dir:          Ref):
+
+    with redirect_logging('fieldmap_pre_topup', outdir=logs_dir),\
+         tempdir(tmp_dir):
+
+        AP_tmp = tmp_dir + '/AP_tmp_'
+        PA_tmp = tmp_dir + '/PA_tmp_'
+
+        choose_best_B0(ctx, AP, AP_bval, total_B0_AP, AP_indices, AP_tmp,
+                       AP_best_index, B0_AP, AP_fieldmap_flag, tmp_dir)
+        choose_best_B0(ctx, PA, PA_bval, total_B0_PA, PA_indices, PA_tmp,
+                       PA_best_index, B0_PA, PA_fieldmap_flag, tmp_dir)
+
+        wrappers.fslmerge("t", B0_AP_PA, B0_AP, B0_PA)
+
+        generate_acqparams(AP, PA, AP_json, PA_json, 1, 1, acqparams)
+
+        # Improvement by Kevin Anderson on 20/12/2017
+        # This guarantees that the number of slices is even
+        im_AP_PA = nib.load(B0_AP_PA)
+        Zslices=im_AP_PA.header['dim'][3]
+
+        if (Zslices % 2) != 0:
+            wrappers.fslroi(B0_AP_PA, B0_AP_PA, 0, -1, 0, -1, 0, Zslices -1)

bip/pipelines/dMRI_fieldmap/fieldmap_topup.py

+#!/usr/bin/env python
+#
+# fieldmap_topup.py - Sub-pipeline with the fieldmap generation with topup
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+# pylint: disable=W0613
+#
+
+import logging
+from fsl import wrappers
+from bip.utils import redirect_logging
+from pipe_tree import In, Out, Ref
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        B0_AP_PA:            In,
+        acqparams:           In,
+        logs_dir:            Ref,
+        fieldmap_out_prefix: Ref,
+        fieldmap_jacout:     Ref,
+        fieldmap_fout:       Out):
+
+    with redirect_logging('fieldmap_topup', outdir=logs_dir):
+        wrappers.topup(imain=B0_AP_PA,
+                       datain=acqparams,
+                       config="b02b0.cnf",
+                       out=fieldmap_out_prefix,
+                       fout=fieldmap_fout,
+                       jacout=fieldmap_jacout,
+                       verbose=True)

bip/pipelines/fMRI_rest/fMRI_rest.py

+#!/usr/bin/env python
+#
+# fMRI_rest.py - Pipeline with the resting state fMRI processing.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import logging
+from bip.utils import redirect_logging
+from bip.pipelines.fMRI_rest import rfMRI_prepare
+from bip.pipelines.fMRI_rest import rfMRI_melodic
+from bip.pipelines.fMRI_rest import rfMRI_fix
+from bip.pipelines.fMRI_rest import rfMRI_netmats
+
+log = logging.getLogger(__name__)
+
+def add_to_pipeline(ctx, pipe, tree, targets):
+
+    logs_dir=tree.get('logs_dir')
+
+    with redirect_logging('pipe_fMRI_task', outdir=logs_dir):
+        pipe(rfMRI_prepare.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        targets.append('rfMRI_fsf')
+        pipe(rfMRI_melodic.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        targets.append('rfMRI_ica')
+        pipe(rfMRI_fix.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        targets.append('filtered_func_data_clean_stdmask')
+        pipe(rfMRI_netmats.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        targets.append('node_amplitudes')
+
+    return pipe, targets

bip/pipelines/fMRI_rest/rfMRI_fix.py

+#!/usr/bin/env python
+#
+# rfMRI_fix.py - Sub-pipeline with the FIX processing of the rest fMRI.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import os
+import logging
+from shutil import copyfile
+from pipe_tree import In, Out, Ref
+from fsl import wrappers
+from bip.utils import redirect_logging
+from pyfix import extract, legacy, io, clean, classify
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        T1_fast_pveseg: In,
+        rfMRI_ica:                        In,
+        rfMRI_example_func2standard_warp: In,
+        rfMRI_example_func:               In,
+        rfMRI_standard:                   In,
+        filtered_func_data:               In,
+        melodic_mix:                      In,
+        rfMRI_mc_parameters:              In,
+        logs_dir:                         Ref,
+        rfMRI_reg_standard_dir:           Ref,
+        rfMRI_highres_pveseg:             Out,
+        rfMRI_standard2example_func_warp: Out,
+        fix_dir:                          Out,
+        fix_features:                     Out,
+        filtered_func_data_clean:         Out,
+        filtered_func_data_clean_std:     Out,
+        filtered_func_data_clean_stdmask: Out,
+        fix4melview:                      Out):
+
+    with redirect_logging('rfMRI_fix', outdir=logs_dir):
+
+        # Generate files needed to extract features
+        if not os.path.exists(rfMRI_highres_pveseg):
+            copyfile(T1_fast_pveseg, rfMRI_highres_pveseg)
+        if not os.path.exists(rfMRI_standard2example_func_warp):
+            wrappers.invwarp(ref=rfMRI_example_func,
+                             warp=rfMRI_example_func2standard_warp,
+                             out=rfMRI_standard2example_func_warp)
+
+        # FIX: Extract features for subject
+        if not os.path.exists(fix_features):
+            d = extract.FixData.from_melodic_dir(rfMRI_ica, fixdir = fix_dir)
+
+            f = extract.extract_features(data=d,
+                                      features=legacy.LEGACY_FEATURE_EXTRACTORS)
+            io.save_features(os.path.join(fix_dir, 'features'), f)
+        else:
+            f = io.read_features(os.path.join(fix_dir, 'features'))
+
+        # FIX: Classify components for subject
+        clf = io.load_model(ctx.get_data('/fix/UKBB_model_2.pyfix_model'))
+        y_pred_pb = classify.classify(f, clf)[1]
+        y_pred = legacy.prob_threshold(y_pred_pb, threshold=0.2)
+
+        fin_list=[]
+        with open(fix4melview, 'wt', encoding="utf-8") as f:
+            for ind, pred in enumerate(y_pred):
+                # Opposite of the convention in MelView
+                if not pred:
+                    f.write(str(ind + 1) + ", Unclassified Noise, True\n")
+                    fin_list.append(ind + 1)
+                else:
+                    f.write(str(ind + 1) + ", Signal, False\n")
+            f.write(f'{fin_list}\n')
+
+        # FIX: Clean
+        clean.apply(func_fname=filtered_func_data,
+                    icamix=melodic_mix,
+                    noise_idx=y_pred,
+                    func_clean_fname=filtered_func_data_clean,
+                    do_motion_regression=True,
+                    motion_parameters=rfMRI_mc_parameters)
+
+        # Create this directory in case it does not exist. This may
+        # happen due to the overwrite function issue in Melodic
+        if not os.path.exists(rfMRI_reg_standard_dir):
+            os.mkdir(rfMRI_reg_standard_dir)
+
+        # Generate the data in standard space
+        wrappers.applywarp(ref=rfMRI_standard, src=filtered_func_data_clean,
+                           out=filtered_func_data_clean_std,
+                           warp=rfMRI_example_func2standard_warp,
+                           interp="spline")
+        mask = ctx.get_data('MNI/MNI152_T1_2mm_brain_mask_bin.nii.gz')
+        wrappers.fslmaths(filtered_func_data_clean_std).mas(mask).run(filtered_func_data_clean_std)
+        wrappers.fslmaths(filtered_func_data_clean_std).Tstd().bin().run(filtered_func_data_clean_stdmask)

bip/pipelines/fMRI_rest/rfMRI_melodic.py

+#!/usr/bin/env python
+#
+# rfMRI_melodic.py - Sub-pipeline with FSL's MELODIC tool
+#                    applied to the rest fMRI.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+# pylint: disable=W0613
+#
+
+import os
+import shutil
+import logging
+from fsl import wrappers
+from bip.utils import redirect_logging
+from pipe_tree import In, Out, Ref
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        rfMRI_fsf:                        In,
+        logs_dir:                         Ref,
+        rfMRI_ica:                        Out,
+        rfMRI_example_func2standard_warp: Out,
+        rfMRI_example_func:               Out,
+        rfMRI_standard:                   Out,
+        filtered_func_data:               Out,
+        melodic_mix:                      Out,
+        rfMRI_mc_parameters:              Out):
+
+    with redirect_logging('rfMRI_melodic', outdir=logs_dir):
+        wrappers.feat(fsf=rfMRI_fsf)
+
+    # The ending of this script is needed because the overwrite output
+    # function in MELODIC is not working, and therefore, the output is
+    # generated in a folder called "rfMI+.ica". Once this function works,
+    # we can remove these last 8 lines
+    name_dir_alt = rfMRI_ica.split(".")
+    name_dir_alt[-2] = name_dir_alt[-2] + "+"
+    name_dir_alt = ".".join(name_dir_alt)
+
+    if os.path.exists(name_dir_alt):
+        if os.path.exists(rfMRI_ica):
+            shutil.rmtree(rfMRI_ica)
+        os.rename(name_dir_alt, rfMRI_ica)

bip/pipelines/fMRI_rest/rfMRI_netmats.py

+#!/usr/bin/env python
+#
+# rfMRI_netmats.py - Sub-pipeline with the netmats processing of the rest fMRI.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import os
+import logging
+from fsl import wrappers
+from pipe_tree import In, Out, Ref, Var
+from bip.utils import redirect_logging
+from bip.pipelines.fMRI_rest.rfMRI_netmats_fnc import ICA_dr
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        netmats_dim:                  Var,
+        filtered_func_data:           In,
+        filtered_func_data_clean_std: In,
+        logs_dir:                     Ref,
+        netmats_dir:                  Out,
+        dr_stage1:                    Out,
+        full_corr:                    Out,
+        partial_corr:                 Out,
+        node_amplitudes:              Out):
+
+    D = netmats_dim.value
+
+    with redirect_logging('rfMRI_netmats_d' + D, outdir=logs_dir):
+
+        MNI_2mm_mask = ctx.get_standard("MNI152_T1_2mm_brain_mask.nii.gz")
+        group_IC = ctx.get_data("/netmats/melodic_IC_" + D + ".nii.gz")
+
+
+        # Create this directory in case it does not exist. This may
+        # happen due to the overwrite function issue in Melodic
+        if not os.path.exists(netmats_dir):
+            os.mkdir(netmats_dir)
+
+
+        wrappers.fsl_glm(input=filtered_func_data_clean_std,
+                         design=group_IC, out=dr_stage1, demean=True,
+                         mask = MNI_2mm_mask)
+
+        # Implement in python
+        ICA_dr(filtered_func_data, dr_stage1, D, full_corr, partial_corr,
+               node_amplitudes)

bip/pipelines/fMRI_rest/rfMRI_netmats_fnc.py

+#!/usr/bin/env python
+#
+# rfMRI_netmats_fnc.py - Functions for the netmats processing of the rest fMRI.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+# pylint: disable=E1130
+#
+
+import numpy as np
+import nibabel as nib
+
+def nets_load(timeseries_name,tr,varnorm,numvols=0):
+
+    timeseries = np.loadtxt(timeseries_name)
+    gn = timeseries.shape[0]
+    gd = timeseries.shape[1]
+
+    ts={}
+
+    ts['NtimepointsPerSubject'] = numvols
+
+    if ts['NtimepointsPerSubject'] == 0:
+        ts["NtimepointsPerSubject"] = gn
+        TS=np.full(timeseries.shape, np.nan)
+
+    TS = timeseries
+    TS = TS - np.nanmean(TS)
+
+    if varnorm == 1:
+        TS = TS / np.nanstd(TS)
+
+    ts["ts"]          = TS
+    ts["tr"]          = tr
+    ts["Nsubjects"]   = 1
+    ts["Nnodes"]      = gd
+    ts["NnodesOrig"]  = ts["Nnodes"]
+    ts["Ntimepoints"] = gn
+    ts["DD"]          = [*range(ts["Nnodes"])]
+    ts["UNK"]         = []
+
+    return ts
+
+def nets_tsclean(ts,aggressive):
+    ts["NnodesOrig"]=ts["Nnodes"]
+
+    nongood = [x for x in range(ts["Nnodes"]) if x not in ts["DD"]]
+    bad = [x for x in nongood if x not in ts["UNK"]] # Only bad components
+
+    # Good components
+    goodTS=ts["ts"][:,ts["DD"]]
+
+    # Bad components
+    badTS=ts["ts"][:,bad]
+
+    if aggressive == 1:
+        ts["ts"] = goodTS - np.dot(badTS, np.dot(np.linalg.pinv(badTS),goodTS))
+    else:
+        ts["ts"] = goodTS
+
+    ts["Nnodes"] = ts["ts"].shape[1]
+
+    return ts
+
+#netmats2 = nets_netmats(ts, -r2zPARTIAL, 'ridgep', 0.5)
+def nets_netmats(ts, do_rtoz, method, arg_method=None):
+
+    N=ts["Nnodes"]
+    netmats=[]
+    grot=ts["ts"]
+
+    # We are keeping just the amplitudes
+    if method.lower() in ['corr','correlation']:
+        grot=np.corrcoef(grot,rowvar=False)
+        np.fill_diagonal(grot, 0)
+
+    elif method.lower() in ['ridgep']:
+        grot=np.cov(grot, rowvar=False)
+        grot=grot/np.sqrt(np.nanmean(np.diag(grot) * np.diag(grot)))
+        if not arg_method:
+            rho=0.1
+        else:
+            rho=arg_method
+            if rho<0:
+                rho = -(rho)
+
+        grot=-(np.linalg.inv(grot+rho*np.eye(N)))
+        val0 = np.sqrt(np.abs(np.diag(grot)))
+        val2 = np.tile(val0, (N,1))
+        val1 = np.transpose(val2)
+        grot = (grot / val1) / val2
+
+        np.fill_diagonal(grot, 0)
+
+    netmats=grot.reshape(1,N*N, order='F').copy()
+    netmats=0.5*np.log((1+netmats)/(1-netmats))* (-do_rtoz)
+
+    return netmats
+
+def ICA_dr(func, timeseries, D, full_corr, partial_corr, node_amplitudes):
+    im_func = nib.load(func)
+    tr      = im_func.header['pixdim'][4]
+    numvols = im_func.header['dim'][4]
+
+    ts = nets_load(timeseries, tr, 0, numvols)
+
+    if D == "25":
+        good_components = [3, 22,23,24]
+        ts["DD"] = [x for x in ts["DD"] if x not in good_components]
+        r2zFULL=10.6484
+        r2zPARTIAL=10.6707
+    else:
+        good_components =  [0, 43, 46, 50, 53, 54, 55, 58, 60, 61]
+        good_components += [*range(64,92)] + [*range(93,100)]
+        ts["DD"] = [x for x in ts["DD"] if x not in good_components]
+        r2zFULL    = 19.7177
+        r2zPARTIAL = 18.8310
+
+    ts = nets_tsclean(ts,1)
+
+    netmats1 = nets_netmats(ts, -r2zFULL,    'corr')
+    netmats2 = nets_netmats(ts, -r2zPARTIAL, 'ridgep', 0.5)
+
+    tmp = np.reshape(netmats1,(ts["Nnodes"],ts["Nnodes"]), order="F")
+    data = tmp[np.where(np.triu(np.ones([ts["Nnodes"],ts["Nnodes"]]),k=1) ==1)]
+    with open(full_corr, 'wt', encoding="utf-8") as f:
+        np.savetxt(f, data, fmt="%14.8f", newline=" ")
+
+    tmp = np.reshape(netmats2,(ts["Nnodes"],ts["Nnodes"]), order="F")
+    data = tmp[np.where(np.triu(np.ones([ts["Nnodes"],ts["Nnodes"]]),k=1) ==1)]
+    with open(partial_corr, 'wt', encoding="utf-8") as f:
+        np.savetxt(f, data, fmt="%14.8f", newline=" ")
+
+    data = np.nanstd(ts["ts"], axis=0)
+    with open(node_amplitudes, 'wt', encoding="utf-8") as f:
+        np.savetxt(f, data, fmt="%14.8f", newline=" ")