diff --git a/docs/Pipeline_Branches.md b/docs/Pipeline_Branches.md new file mode 100644 index 0000000..4e3cfa5 --- /dev/null +++ b/docs/Pipeline_Branches.md @@ -0,0 +1,49 @@ +There are many tasks in our Luigi pipeline defined in [different modules](../workflows): + +* struct_pipe.py +* dwi_pipe.py +* fs2dwi_pipe.py + +Some of them can be run standalone. Others require a task to be run beforehand. Such various routes are outlined below: + +--- + +#### T1 + T2 templates + + StructMask, Freesurfer + +#### dwi + T2 templates + + EddyEpi, Ukf, Wma800 + +Provided that `Freesurfer` task was run: + + EddyEpi, Ukf, Fs2Dwi, Wmql, Wmqlqc, TractMeasures + +--- + +#### Only T1 template + + StructMask, Freesurfer + +#### dwi+T1 template + + SynB0, Ukf, Wma800 + +#### ap_pa_template + + TopupEddy, Ukf, Wma800 + +--- + +#### Only dwi template + + DwiAlign, GibbsUn, CnnMask, PnlEddy, FslEddy, Ukf, Wma800 + +Provided that `Freesurfer` task was run: + + DwiAlign, GibbsUn, CnnMask, PnlEddy, FslEddy, Ukf, Fs2Dwi, Wmql, Wmqlqc, TractMeasures + +Provided that `GibbsUn` task was run: + + HcpPipe, Ukf, Wma800 diff --git a/docs/SynB0-Wma800.png b/docs/SynB0-Wma800.png new file mode 100644 index 0000000..d445f2f Binary files /dev/null and b/docs/SynB0-Wma800.png differ diff --git a/docs/TUTORIAL.md b/docs/TUTORIAL.md index 34bd5b9..47bfe51 100644 --- a/docs/TUTORIAL.md +++ b/docs/TUTORIAL.md @@ -23,7 +23,9 @@ Table of Contents * [Correct eddy and epi](#correct-eddy-and-epi) * [EddyEpi](#eddyepi) * [FslEddy and PnlEddy](#fsleddy-and-pnleddy) + * [SynB0](#synb0) * [TopupEddy](#topupeddy) + * [HcpPipe](#hcppipe) * [UKFTractography](#ukftractography) * [Wma800](#wma800) * [Fs2Dwi pipeline](#fs2dwi-pipeline) @@ -384,7 +386,9 @@ Different methods for eddy and epi corrections are described below: * PnlEddy + PnlEpi (EddyEpi) * FslEddy + PnlEpi (EddyEpi) +* SynB0 * TopupEddy +* HcpPipe ### EddyEpi @@ -459,6 +463,47 @@ The mask and baseline image provided to `FslEddy` are approximated as eddy corre be fed into later task FslEddyEpi. +### SynB0 + +When an axial T2w image is absent and there is only AP or PA acquisition, we make use of [Synb0-DISCO](https://github.com/MASILab/Synb0-DISCO) +program to perform EPI distortion correction. The Synb0-DISCO program is a replacement of topup. It creates a distortion corrected B0 from a single +AP or PA acquisition. It also provides the usual topup correction parameters that are used to perform eddy correction later. +However, Synb0-DISCO requires a T1w image. + +![](https://github.com/pnlbwh/luigi-pnlpipe/blob/0d2229e1f27fcc6b02825e9a147c0146924eae43/docs/SynB0-Wma800.png) + + +The configuration file used for `SynB0` task is same as the above [dwi_pipe_params.cfg](../params/dwi_pipe_params.cfg). Specifically, `SynB0` requires +the following parameters definition: + +```cfg +## [SynB0] ## +acqp: /path/to/acqp.txt +index: /path/to/index.txt +``` + +Where `acqp.txt` just: + +``` +0 1 0 0.05 +0 1 0 0 +``` + +Run it: + +```bash +export LUIGI_CONFIG_PATH=/path/to/dwi_pipe_params.cfg + +exec/ExecuteTask --task SynB0 \ +--bids-data-dir /data/pnl/Collaborators/EDCRP/1110/1110_ses-002_T1_nii/rawdata/ \ +-c ne00181 -s 002 \ +--dwi-template sub-*/ses-*/dwi/*_dwi.nii.gz \ +--t1-template sub-*/ses-*/anat/*_T1w.nii.gz +``` + +Note that, `SynB0` task is run in coordination with [_synb0_eddy.sh](_synb0_eddy.sh) . + + ### TopupEddy When two opposing acquisitions--AP and PA are available such as in Human Connectom Project (HCP), @@ -506,6 +551,22 @@ for corrected data is determined as follows: *_dir-214_* +### HcpPipe + +This task supersedes *TopupEddy* task. The *HcpPipe* way of processing opposing pair of DWIs is inherited +from [Washington-University/HCPpipelines](https://github.com/Washington-University/HCPpipelines). +PNL introduced some improvements in that pipeline including: + + * Use [CNN brain masking](https://github.com/pnlbwh/CNN-Diffusion-MRIBrain-Segmentation) program to generate topup mask + * Replace outliers (`--repol` flag) for >500 b-shells + +Since external processing is involved in this task, it is run via a [separate script](https://github.com/pnlbwh/luigi-pnlpipe/blob/b311000073047b40db9b7b8c1752cdbba883aa6c/workflows/hcp_pnl_topup.lsf) in three steps: + + * Gibbs unringing of opposing pair of DWIs via Luigi + * HCP pipeline via shell scripts + * Creation of soft links in `sub-*/ses-*/dwi/` directory according to BIDS convention via Luigi + +This *HcpPipe* task has been used to process *HCP-EP* data. ## UKFTractography diff --git a/workflows/ExecuteTask.py b/workflows/ExecuteTask.py index ffc933c..cc546fe 100755 --- a/workflows/ExecuteTask.py +++ b/workflows/ExecuteTask.py @@ -6,7 +6,7 @@ from _define_outputs import IO from struct_pipe import StructMask, Freesurfer from dwi_pipe import DwiAlign, GibbsUn, CnnMask, \ - PnlEddy, FslEddy, TopupEddy, HcpPipe, EddyEpi, Ukf + PnlEddy, FslEddy, SynB0, TopupEddy, HcpPipe, EddyEpi, Ukf, Wma800 from fs2dwi_pipe import Fs2Dwi, Wmql, Wmqlqc, TractMeasures from scripts.util import abspath, isfile, pjoin, LIBDIR from os import getenv, stat, remove @@ -56,9 +56,9 @@ def _rm_tempfiles(names): default= argparse.SUPPRESS, choices=['StructMask', 'Freesurfer', 'DwiAlign', 'GibbsUn', 'CnnMask', - 'PnlEddy', 'FslEddy', 'TopupEddy', 'HcpPipe', + 'PnlEddy', 'FslEddy', 'SynB0', 'TopupEddy', 'HcpPipe', 'EddyEpi', - 'Ukf', + 'Ukf', 'Wma800', 'Fs2Dwi', 'Wmql', 'Wmqlqc', 'TractMeasures']) parser.add_argument('--num-workers', type=int, default=1, help='number of Luigi workers') @@ -93,158 +93,125 @@ def _rm_tempfiles(names): if args.t2_template: if args.task=='StructMask': - jobs.append(StructMask(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - struct_template=args.t2_template)) + jobs.append(StructMask( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + struct_template=args.t2_template)) elif args.task=='Freesurfer': - jobs.append(Freesurfer(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - t1_template=args.t1_template, - t2_template=args.t2_template)) + jobs.append(Freesurfer( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + t1_template=args.t1_template, + t2_template=args.t2_template)) - elif args.task=='EddyEpi': - jobs.append(eval(args.task)(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template, - struct_template=args.t2_template)) - - - - # Ukf task does not have pa_ap_template because - # when axt2 is available, pa_ap acquisition should be unavailable - # in other words, PnlEpi and TopupEddy are mutually exclusive - elif args.task=='Ukf': - jobs.append(Ukf(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template, - struct_template=args.t2_template)) - - - elif args.task=='Fs2Dwi': - jobs.append(Fs2Dwi(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template, - struct_template=args.t2_template)) - - elif args.task=='Wmql': - jobs.append(Wmql(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template, - struct_template=args.t2_template)) - - elif args.task=='TractMeasures': - jobs.append(TractMeasures(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template, - struct_template=args.t2_template)) - - - - elif args.task=='Wmqlqc': - jobs.append(Wmqlqc(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template, - struct_template=args.t2_template)) - + jobs.append(EddyEpi( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + dwi_template=args.dwi_template, + struct_template=args.t2_template)) + + + # This Ukf task does not have pa_ap_template because + # when AXT2 is available, pa_ap acquisition should be unavailable. + # In other words, PnlEpi and TopupEddy are mutually exclusive + elif args.task in ['Ukf','Wma800']: + jobs.append(eval(args.task)( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + dwi_template=args.dwi_template, + struct_template=args.t2_template)) + + elif args.task in ['Fs2Dwi','Wmql','TractMeasures','Wmqlqc']: + jobs.append(eval(args.task)( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + dwi_template=args.dwi_template, + struct_template=args.t2_template)) + + # just t1_template else: if args.task=='StructMask': - jobs.append(StructMask(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - struct_template=args.t1_template)) + jobs.append(StructMask( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + struct_template=args.t1_template)) elif args.task=='Freesurfer': - jobs.append(Freesurfer(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - t1_template=args.t1_template)) - - + jobs.append(Freesurfer( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + t1_template=args.t1_template)) elif args.task in ['DwiAlign','GibbsUn','CnnMask','PnlEddy','FslEddy','HcpPipe']: - jobs.append(eval(args.task)(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template)) - + jobs.append(eval(args.task)( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + dwi_template=args.dwi_template)) + + elif args.task=='SynB0': + jobs.append(SynB0( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + dwi_template=args.dwi_template, + struct_template=args.t1_template)) elif args.task=='TopupEddy': - jobs.append(TopupEddy(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - pa_ap_template=args.dwi_template)) - - # Ukf task has both dwi_template and pa_ap_template - # because a user may want to run {PnlEddy,FslEddy} or TopupEddy - elif args.task=='Ukf': - jobs.append(Ukf(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template, - pa_ap_template=args.dwi_template)) - - - elif args.task=='Fs2Dwi': - jobs.append(Fs2Dwi(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template)) - - - elif args.task == 'Wmql': - jobs.append(Wmql(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template)) - - - elif args.task=='Wmqlqc': - jobs.append(Wmqlqc(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template)) - - - elif args.task=='TractMeasures': - jobs.append(TractMeasures(bids_data_dir=args.bids_data_dir, - derivatives_dir=derivatives_dir, - id=id, - ses=ses, - dwi_template=args.dwi_template)) + jobs.append(TopupEddy( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + pa_ap_template=args.dwi_template)) + + + # This Ukf task has both dwi_template and pa_ap_template + # because a user may want to run {PnlEddy,FslEddy,HcpPipe} or TopupEddy + elif args.task in ['Ukf','Wma800']: + jobs.append(eval(args.task)( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + dwi_template=args.dwi_template, + pa_ap_template=args.dwi_template, + struct_template=args.t1_template)) + + elif args.task in ['Fs2Dwi','Wmql','TractMeasures','Wmqlqc']: + jobs.append(eval(args.task)( + bids_data_dir=args.bids_data_dir, + derivatives_dir=derivatives_dir, + id=id, + ses=ses, + dwi_template=args.dwi_template)) + build(jobs, workers=args.num_workers) - print('Removing temporary provenance files') + # print('Removing temporary provenance files') _rm_tempfiles(glob(pjoin(gettempdir(), 'hashes-*.txt'))) _rm_tempfiles(glob(pjoin(gettempdir(), 'env-*.yml'))) diff --git a/workflows/_synb0_eddy.sh b/workflows/_synb0_eddy.sh new file mode 100755 index 0000000..bac2046 --- /dev/null +++ b/workflows/_synb0_eddy.sh @@ -0,0 +1,90 @@ +#!/bin/bash + +# for time profiling +date + +unring_prefix=${1//.nii.gz/} +b0=$2 +T1=$3 +eddy_prefix=${4//.nii.gz/} +eddy_mask=$5 +eddy_bse=$6 +ACQPARAMS=$7 +INDEX=$8 + + +pushd . +SES_FOLDER=$(dirname $(dirname $unring_prefix)) +cd $SES_FOLDER +mkdir -p INPUTS OUTPUTS + +cp $b0 INPUTS/b0.nii.gz +cp $T1 INPUTS/T1.nii.gz + +cp $ACQPARAMS INPUTS/ +cp $INDEX INPUTS/ + + +echo "3. run synb0 container" +TMPDIR=$HOME/tmp/ +mkdir -p $TMPDIR +if [ ! -f OUTPUTS/b0_all_topup.nii.gz ] +then + TMPDIR=$TMPDIR \ + singularity run -B INPUTS/:/INPUTS -B OUTPUTS/:/OUTPUTS \ + -B ${NEW_SOFT_DIR}/fs7.1.0/license.txt:/extra/freesurfer/license.txt \ + ${NEW_SOFT_DIR}/containers/synb0-disco_v3.0.sif --stripped +fi + +echo "4. create mask of topup (synb0) corrected b0" +# CNN method +_caselist=$(mktemp --suffix=.txt) +realpath OUTPUTS/b0_all_topup.nii.gz > $_caselist +echo "0 0" > OUTPUTS/b0_all_topup.bval +dwi_masking.py -i $_caselist -f ${NEW_SOFT_DIR}/CNN-Diffusion-MRIBrain-Segmentation/model_folder +mask=`ls OUTPUTS/*-multi_BrainMask.nii.gz` +rm $_caselist + +if [ -z $mask ] +then + echo topup mask creation failed + exit 1 +fi + + +_eddy_out=OUTPUTS/$(basename $eddy_prefix) +echo "5. run eddy_cuda" +eddy_cuda \ + --imain=${unring_prefix}.nii.gz \ + --bvecs=${unring_prefix}.bvec \ + --bvals=${unring_prefix}.bval \ + --mask=$mask \ + --topup=OUTPUTS/topup \ + --acqp=INPUTS/acqparams.txt \ + --index=INPUTS/index.txt \ + --repol --data_is_shelled --verbose \ + --out=${_eddy_out} || { exit 1; } + + + +echo "6. organize outputs" +# provide masked eddy_out for clarity, quality, and convenience +fslmaths ${_eddy_out}.nii.gz -mul $mask ${_eddy_out}.nii.gz + +mv ${_eddy_out}.nii.gz ${eddy_prefix}.nii.gz +mv ${_eddy_out}.eddy_rotated_bvecs ${eddy_prefix}.bvec +cp ${unring_prefix}.bval ${eddy_prefix}.bval + +mv $mask $eddy_mask +mv OUTPUTS/b0_all_topup_bse.nii.gz $eddy_bse + + + +echo "Luigi-SynB0-Eddy pipeline has completed" +echo "See outputs at $PWD/dwi/" + +popd + +# for time profiling +date + diff --git a/workflows/dwi_pipe.py b/workflows/dwi_pipe.py index d4b85a6..6b5ec25 100755 --- a/workflows/dwi_pipe.py +++ b/workflows/dwi_pipe.py @@ -12,7 +12,7 @@ from time import sleep import re -from struct_pipe import StructMask +from struct_pipe import StructMask, N4BiasCorrect from _task_util import _mask_name from scripts.util import N_PROC, B0_THRESHOLD, BET_THRESHOLD, QC_POLL, LIBDIR, \ @@ -87,7 +87,7 @@ def output(self): dwi = local.path(re.sub('_desc-Xc_', '_desc-XcUn_', self.input()['dwi'])) bval = dwi.with_suffix('.bval', depth=2) bvec = dwi.with_suffix('.bvec', depth=2) - + return dict(dwi=dwi, bval=bval, bvec=bvec) @@ -323,6 +323,69 @@ def output(self): +@requires(GibbsUn,CnnMask,N4BiasCorrect) +class SynB0(Task): + + mask_qc= BoolParameter(default=False) + acqp = Parameter() + index = Parameter() + + + def run(self): + + # synb0 wrapper + DIR= abspath(dirname(__file__)) + cmd = (' ').join([f'{DIR}/_synb0_eddy.sh', + self.input()[0]['dwi']._path, + self.input()[1]['bse']._path, + self.input()[2]['n4corr']._path, + self.output()['dwi']._path, + self.output()['mask'], + self.output()['bse'], + self.acqp, + self.index]) + p = Popen(cmd, shell=True) + p.wait() + + version_file= self.output()['dwi'].dirname.join('fsl_version.txt') + check_call(f'eddy_openmp 2>&1 | grep Part > {version_file}', shell= True) + + + write_provenance(self, self.output()['dwi']) + + # self.dwi= self.output()['dwi'] + # yield self.clone(BseExtract) + + + + def output(self): + + eddy_epi_prefix= self.input()[0]['dwi'].rsplit('_dwi.nii.gz')[0] + eddy_epi_prefix= eddy_epi_prefix.replace('_acq-PA','') + eddy_epi_prefix= eddy_epi_prefix.replace('_acq-AP','') + eddy_epi_prefix+= 'EdEp' + + dwi = local.path(eddy_epi_prefix+ '_dwi.nii.gz') + bval = dwi.with_suffix('.bval', depth= 2) + bvec = dwi.with_suffix('.bvec', depth= 2) + + + mask_prefix = dwi.rsplit('_desc-')[0] + desc= re.search('_desc-(.+?)_dwi.nii.gz', basename(dwi)).group(1) + desc= 'dwi'+ desc + mask_prefix= mask_prefix+ '_desc-'+ desc + mask = local.path(mask_prefix+ '_mask.nii.gz') + + bse_prefix= dwi._path + desc= re.search('_desc-(.+?)_dwi.nii.gz', bse_prefix).group(1) + desc= 'dwi'+ desc + bse= local.path(bse_prefix.split('_desc-')[0]+ '_desc-'+ desc+ '_bse.nii.gz') + + + return dict(dwi=dwi, bval=bval, bvec=bvec, bse=bse, mask=mask) + + + @inherits(FslEddy, PnlEddy, StructMask, BseExtract) class EddyEpi(Task): debug = BoolParameter(default=False) @@ -405,9 +468,9 @@ class TopupEddy(Task): config = Parameter(default=pjoin(LIBDIR, 'scripts', 'eddy_config.txt')) useGpu = BoolParameter(default=False) - numb0 = Parameter(default=1) + numb0 = Parameter(default='1') whichVol = Parameter(default='1') - scale = Parameter(default=2) + scale = Parameter(default='2') TopupOutDir= Parameter(default='fsl_eddy') @@ -595,7 +658,7 @@ def output(self): -@inherits(PnlEddy, FslEddy, EddyEpi, TopupEddy, HcpPipe) +@inherits(PnlEddy, FslEddy, SynB0, EddyEpi, TopupEddy, HcpPipe) class Ukf(Task): ukf_params = Parameter(default='') @@ -609,6 +672,8 @@ def requires(self): return self.clone(PnlEddy) elif self.eddy_epi_task=='fsleddy': return self.clone(FslEddy) + elif self.eddy_epi_task=='synb0': + return self.clone(SynB0) elif self.eddy_epi_task=='eddyepi': return self.clone(EddyEpi) elif self.eddy_epi_task=='topupeddy': @@ -654,7 +719,7 @@ class Wma800(Task): wma_cleanup= IntParameter(default=0) def run(self): - + outDir = self.input().dirname.join('wma800') cmd = (' ').join(['wm_apply_ORG_atlas_to_subject.sh', '-i', self.input(), '-a', self.atlas, @@ -664,12 +729,19 @@ def run(self): f'-n {self.wma_nproc}', f'-c {self.wma_cleanup}', '-d 1', - '-o', self.output()]) + '-o', outDir]) p = Popen(cmd, shell=True) p.wait() - - write_provenance(self) + + write_provenance(self, outDir) def output(self): - return self.input().dirname.join('wma800') + prefix= self.input().dirname.join('wma800',self.input().basename.split('.vtk')[0], + 'FiberClustering/SeparatedClusters') + + clusters=[] + for region in 'commissural left_hemisphere right_hemisphere'.split(): + clusters.append( local.path(f'{prefix}/diffusion_measurements_{region}.csv') ) + + return clusters diff --git a/workflows/synb0_eddy.sh b/workflows/synb0_eddy.sh index 386096a..a4bf611 100755 --- a/workflows/synb0_eddy.sh +++ b/workflows/synb0_eddy.sh @@ -2,6 +2,7 @@ # for time profiling date +START_TIME=$(date +%s) # User will edit only this block ========================================================= caselist=$1 @@ -35,7 +36,6 @@ fi echo "1. run Luigi pipeline and prepare DWI for synb0 container" export LUIGI_CONFIG_PATH -# source /rfanfs/pnl-zorro/software/pnlpipe3/bashrc3-gpu \ /rfanfs/pnl-zorro/software/pnlpipe3/luigi-pnlpipe/exec/ExecuteTask --task CnnMask \ --bids-data-dir $BIDS_DATA_DIR \ --dwi-template "$DWI_TEMPLATE" \ @@ -64,7 +64,6 @@ unring_prefix=${_unring_prefix//.nii.gz} unring_mask=`ls dwi/sub-${c}_ses-${s}_*desc-dwiXcUnCNN_mask.nii.gz` if [ ! -f INPUTS/b0.nii.gz ] then - # source /rfanfs/pnl-zorro/software/pnlpipe3/bashrc3-gpu && \ fslmaths ${unring_prefix}.nii.gz -mul $unring_mask ${unring_prefix}.nii.gz && \ bse.py -i ${unring_prefix}.nii.gz -o INPUTS/b0.nii.gz fi @@ -98,16 +97,9 @@ echo "4. create mask of topup (synb0) corrected b0" _caselist=$(mktemp --suffix=.txt) realpath OUTPUTS/b0_all_topup.nii.gz > $_caselist echo "0 0" > OUTPUTS/b0_all_topup.bval -# source /rfanfs/pnl-zorro/software/pnlpipe3/bashrc3-gpu && \ dwi_masking.py -i $_caselist -f ${NEW_SOFT_DIR}/CNN-Diffusion-MRIBrain-Segmentation/model_folder mask=`ls OUTPUTS/*-multi_BrainMask.nii.gz` rm $_caselist -# BET method -# cd OUTPUTS/ -# fslroi b0_all_topup.nii.gz _b0.nii.gz 0 1 -# bet _b0.nii.gz b0_all_topup -m -n -# mask=`realpath b0_all_topup_mask.nii.gz` -# cd .. if [ -z $mask ] then @@ -123,7 +115,6 @@ eddy_out=OUTPUTS/sub-${c}_ses-${s}_dir-${dir}_desc-XcUnEdEp_dwi # so omit masking at this step # fslmaths ${unring_prefix}.nii.gz -mul $mask ${unring_prefix}.nii.gz echo "5. run eddy_cuda" -# source /rfanfs/pnl-zorro/software/pnlpipe3/bashrc3-gpu-cuda-10.2 && \ eddy_cuda \ --imain=${unring_prefix}.nii.gz \ --bvecs=${unring_prefix}.bvec \ @@ -145,7 +136,11 @@ bids_prefix=dwi/sub-${c}_ses-${s}_dir-${dir}_desc-XcUnEdEp_dwi mv ${eddy_out}.nii.gz ${bids_prefix}.nii.gz mv ${eddy_out}.eddy_rotated_bvecs ${bids_prefix}.bvec cp ${unring_prefix}.bval ${bids_prefix}.bval -mv $mask dwi/sub-${c}_ses-${s}_dir-${dir}_desc-dwiXcUnEdEp_mask.nii.gz + +mask_prefix=${bids_prefix//_dwi/} +mask_prefix=${mask_prefix//XcUn/dwiXcUn} +mv $mask ${mask_prefix}_mask.nii.gz +mv OUTPUTS/b0_all_topup_bse.nii.gz ${mask_prefix}_bse.nii.gz @@ -156,4 +151,5 @@ popd # for time profiling date - +ELAPSED_TIME=$(($(date +%s) - $START_TIME)) +echo "Elapsed time: $(($ELAPSED_TIME/60)) min $(($ELAPSED_TIME%60)) sec"