diff --git a/USAGE b/USAGE
index c3407e1..f92f339 100755
--- a/USAGE
+++ b/USAGE
@@ -37,6 +37,8 @@ DESCRIPTION
 
     --bids=/path/to/[bids]              BIDS folder following the BIDS convention.
 
+    --fs=/path/to/[fs]                  Freesurfer output folder.
+
     --bids_config                       BIDS configuration file in json pre-generated by TractoFlow(Read_BIDS process).
                                         This option must be used only if --bids did not work due to an issue in the BIDS
                                         folder.
@@ -109,9 +111,10 @@ OPTIONAL ARGUMENTS (current value)
 
 --fa                                        Initial FA threshold to compute the frf ($fa).
 --min_fa                                    Minimum FA threshold to compute the frf ($min_fa).
+--min_nvox                                  Minimum number of voxels to compute the frf ($min_nvox).
 --roi_radius                                Region of interest radius to compute the frf ($roi_radius).
 --set_frf                                   Set manually the frf ($set_frf).
---manual_frf                                FRF set manually ($manual_frf).
+--manual_frf                                FRF set manually (--manual_frf "$manual_frf").
 
 --mean_frf                                  Mean the frf of all subjects ($mean_frf).
                                             USE ONLY IF ALL OF SUBJECTS COME FROM THE SAME SCANNER
@@ -208,7 +211,8 @@ cbrain                                      When this profile is used, Nextflow
 
 ABS                                         When this profile is used, TractoFlow-ABS (Atlas Based Segmentation) is used.
                                             This profile must be used for pathological data.
-                                            The aparc+aseg.nii.gz and wmparc.nii.gz must be in the same space than t1.nii.gz
+                                            The aparc+aseg.nii.gz and wmparc.nii.gz must be in the same space than t1.nii.gz.
+                                            In order to use BIDS, either use the root structure or bids and freesurfer structures.
 
 bundling                                    When this profile is used, it will activate custom tracking parameters to improve recobundle results.
                                             Local tracking will be enable with fa seeding mask and tracking mask.
diff --git a/main.nf b/main.nf
old mode 100755
new mode 100644
index a99a552..6999bde
--- a/main.nf
+++ b/main.nf
@@ -3,6 +3,7 @@
 import groovy.json.*
 
 params.input = false
+params.fs = false
 params.bids = false
 params.bids_config = false
 params.help = false
@@ -45,6 +46,7 @@ if(params.help) {
                 "number_of_tissues":"$params.number_of_tissues",
                 "fa":"$params.fa",
                 "min_fa":"$params.min_fa",
+                "min_nvox":"$params.min_nvox",
                 "roi_radius":"$params.roi_radius",
                 "set_frf":"$params.set_frf",
                 "manual_frf":"$params.manual_frf",
@@ -139,20 +141,24 @@ if (params.input && !(params.bids && params.bids_config)){
                        flat: true) {it.parent.name}
 
     labels_for_reg = Channel
-            .fromFilePairs("$root/**/*{aparc+aseg.nii.gz,wmparc.nii.gz}",
-                           size: 2,
-                           maxDepth:1,
-                           flat: true) {it.parent.name}
+        .fromFilePairs("$root/**/*{aparc+aseg.nii.gz,wmparc.nii.gz}",
+                        size: 2,
+                        maxDepth:1,
+                        flat: true) {it.parent.name}
 
-    data
-        .map{[it, params.readout, params.encoding_direction].flatten()}
+    data.map{[it[0..3], "_", it[4], params.readout, params.encoding_direction].flatten()}
         .into{in_data; check_subjects_number}
 
     Channel
-    .fromPath("$root/**/*rev_b0.nii.gz",
-                    maxDepth:1)
-    .map{[it.parent.name, it]}
-    .into{rev_b0; check_rev_b0}
+        .fromPath("$root/**/*rev_b0.nii.gz",
+                        maxDepth:1)
+        .map{[it.parent.name, it]}
+        .tap{rev_b0_for_topup; check_simple_rev_b0}
+        .map{ [it[0]] }
+        .into{sid_rev_b0_included; sid_rev_b0_included_for_eddy_topup}
+
+    Channel.empty().into{sid_rev_dwi_included; sid_rev_dwi_included_for_topup; sid_rev_dwi_excluded; check_rev_number}
+    Channel.empty().into{complex_rev_b0_for_topup; check_complex_rev_b0}
 }
 else if (params.bids || params.bids_config){
     if (!params.bids_config) {
@@ -162,6 +168,11 @@ else if (params.bids || params.bids_config){
             participant_cleaned = workflow.commandLine.substring(start, workflow.commandLine.indexOf("--", start) == -1 ? workflow.commandLine.length() : workflow.commandLine.indexOf("--", start)).replace("=", "").replace("\'", "")
             log.info "Participants: $participant_cleaned"
         }
+        if (params.fs) {
+            Integer start = workflow.commandLine.indexOf("fs") + "fs".length();
+            freesurfer_path = workflow.commandLine.substring(start, workflow.commandLine.indexOf("--", start) == -1 ? workflow.commandLine.length() : workflow.commandLine.indexOf("--", start)).replace("=", "").replace("\'", "").split('-')[0]
+            log.info "Freesurfer path: $freesurfer_path"
+        }
         log.info "Clean_bids: $params.clean_bids"
         log.info ""
 
@@ -183,12 +194,14 @@ else if (params.bids || params.bids_config){
             script:
             participants_flag =\
             params.participants_label ? '--participants_label ' + participant_cleaned : ""
+            fs_flag =\
+            params.fs ? '--fs ' + freesurfer_path : ""
 
             clean_flag = params.clean_bids ? '--clean ' : ''
 
             """
             scil_validate_bids.py $bids_folder tractoflow_bids_struct.json\
-                --readout $params.readout $participants_flag $clean_flag
+                --readout $params.readout $participants_flag $clean_flag $fs_flag -v
             """
         }
     }
@@ -200,7 +213,12 @@ else if (params.bids || params.bids_config){
     }
 
     ch_in_data = Channel.create()
-    ch_rev_b0 = Channel.create()
+    ch_sid_rev_dwi = Channel.create()
+    ch_sid_rev_b0 = Channel.create()
+    ch_complex_rev_b0 = Channel.create()
+    ch_simple_rev_b0 = Channel.create()
+    labels_for_reg = Channel.create()
+
     bids_struct.map{it ->
     jsonSlurper = new JsonSlurper()
         data = jsonSlurper.parseText(it.getText())
@@ -230,26 +248,54 @@ else if (params.bids || params.bids_config){
                     "using --bids."
                 }
             }
-            sub = [sid, file(item.bval), file(item.bvec), file(item.dwi),
+            sub = [sid, file(item.bval), file(item.bvec), file(item.dwi), "_",
                    file(item.t1), item.TotalReadoutTime, item.DWIPhaseEncodingDir[0]]
             ch_in_data.bind(sub)
 
-            if(item.rev_b0) {
-                sub_rev_b0 = [sid, file(item.rev_b0)]
-                ch_rev_b0.bind(sub_rev_b0)}
-        }
+            if(item.rev_topup) {
+                ch_sid_rev_b0.bind([sid])
+                if(item.topup) {
+                  sub_complex_rev_b0 = [sid, file(item.rev_topup), file(item.topup)]
+                  ch_complex_rev_b0.bind(sub_complex_rev_b0)
+                }
+                else{
+                  sub_simple_rev_b0 = [sid, file(item.rev_topup)]
+                  ch_simple_rev_b0.bind(sub_simple_rev_b0)
+                }
+            }
+            else if(item.rev_dwi){
+                ch_rev_in_data = [sid, file(item.rev_bval), file(item.rev_bvec), file(item.rev_dwi), "_rev_",
+                                    file(item.t1), item.TotalReadoutTime, item.DWIPhaseEncodingDir[0]]
+                ch_sid_rev_dwi.bind([sid])
+                ch_in_data.bind(ch_rev_in_data)
+            }
 
+            if(item.wmparc) {
+                sub_labels_for_reg = [sid, file(item.aparc_aseg), file(item.wmparc)]
+                labels_for_reg.bind(sub_labels_for_reg)
+            }
+        }
+        ch_sid_rev_dwi.close()
+        ch_sid_rev_b0.close()
         ch_in_data.close()
-        ch_rev_b0.close()
+        ch_simple_rev_b0.close()
+        ch_complex_rev_b0.close()
+        labels_for_reg.close()
     }
 
+    ch_sid_rev_dwi.into{sid_rev_dwi_included; sid_rev_dwi_included_for_topup; sid_rev_dwi_excluded; check_rev_number}
+    ch_sid_rev_b0.into{sid_rev_b0_included; sid_rev_b0_included_for_eddy_topup}
     ch_in_data.into{in_data; check_subjects_number}
-    ch_rev_b0.into{rev_b0; check_rev_b0}
+
+    ch_simple_rev_b0.into{rev_b0_for_topup; check_simple_rev_b0}
+    ch_complex_rev_b0.into{complex_rev_b0_for_topup; check_complex_rev_b0}
 }
 else {
     error "Error ~ Please use --input, --bids or --bids_config for the input data."
 }
 
+check_subjects_number.map{[it[0]]}.unique().set{unique_subjects_number}
+
 if (params.sh_fitting && !params.sh_fitting_shells){
     error "Error ~ Please set the SH fitting shell to use."
 }
@@ -286,23 +332,37 @@ if (params.run_pft_tracking && workflow.profile.contains("ABS")){
     error "Error ~ PFT tracking cannot be run with Atlas Based Segmentation (ABS) profile"
 }
 
-if (params.bids && workflow.profile.contains("ABS")){
+if (params.bids && workflow.profile.contains("ABS") && !params.fs){
     error "Error ~ --bids parameter cannot be run with Atlas Based Segmentation (ABS) profile"
 }
 
 (dwi, gradients, t1, readout_encoding) = in_data
-    .map{sid, bvals, bvecs, dwi, t1, readout, encoding -> [tuple(sid, dwi),
-                                        tuple(sid, bvals, bvecs),
+    .map{sid, bvals, bvecs, dwi, rev_flag, t1, readout, encoding -> [tuple(sid, dwi, rev_flag),
+                                        tuple(sid, bvals, bvecs, rev_flag),
                                         tuple(sid, t1),
                                         tuple(sid, readout, encoding)]}
     .separate(4)
 
-t1
+t1.unique()
     .into{t1_for_denoise; t1_for_test_denoise}
 
-check_rev_b0.count().into{ rev_b0_counter; number_rev_b0_for_compare }
+check_complex_rev_b0.concat(check_simple_rev_b0).count().into{rev_b0_counter; number_rev_b0_for_compare}
+
+unique_subjects_number.count().into{number_subj_for_null_check; number_subj_for_compare}
 
-check_subjects_number.count().into{ number_subj_for_null_check; number_subj_for_compare }
+check_rev_number.count().into{number_rev_dwi; rev_dwi_counter}
+
+if (params.eddy_cmd == "eddy_openmp"){
+number_rev_dwi
+    .subscribe{a -> if (a>0)
+    error "Error ~ You have some subjects with a reverse encoding DWI. You MUST add -profile use_cuda with a GPU environnement to be able to analyse these data."}
+}
+
+if (!params.run_topup || !params.run_eddy){
+number_rev_dwi
+    .subscribe{a -> if (a>0)
+    error "Error ~ You have some subjects with a reverse encoding DWI. You MUST run topup and eddy with this kind of acquisition."}
+}
 
 number_subj_for_null_check
 .subscribe{a -> if (a == 0)
@@ -338,17 +398,15 @@ else{
 }
 
 gradients
-    .into{gradients_for_prelim_bet; gradients_for_eddy; gradients_for_topup;
+    .into{gradients_for_prelim_bet; gradients_for_eddy;
+          gradients_for_prepare_topup;
+          gradients_for_prepare_dwi_for_eddy;
           gradients_for_eddy_topup; gradients_for_test_eddy_topup}
 
 readout_encoding
     .into{readout_encoding_for_topup; readout_encoding_for_eddy;
           readout_encoding_for_eddy_topup}
 
-dwi_for_prelim_bet
-    .join(gradients_for_prelim_bet)
-    .set{dwi_gradient_for_prelim_bet}
-
 process README {
     cpus 1
     publishDir = params.Readme_Publish_Dir
@@ -373,12 +431,19 @@ process README {
     """
 }
 
+dwi_for_prelim_bet
+    .map{ [it[0] + it[2]] + it }
+    .join(gradients_for_prelim_bet.map{ [it[0] + it[3], it[1], it[2]] })
+    .map{ it[1..-1] }
+    .into{dwi_gradient_for_prelim_bet;toto}
+
 process Bet_Prelim_DWI {
     cpus 2
 
     input:
-    set sid, file(dwi), file(bval), file(bvec) from dwi_gradient_for_prelim_bet
+    set sid, file(dwi), val(rev), file(bval), file(bvec) from dwi_gradient_for_prelim_bet
     val(rev_b0_count) from rev_b0_counter
+    val(rev_dwi_count) from rev_dwi_counter
 
     output:
     set sid, "${sid}__b0_bet_mask_dilated.nii.gz" into\
@@ -387,7 +452,7 @@ process Bet_Prelim_DWI {
     file "${sid}__b0_bet_mask.nii.gz"
 
     when:
-    params.run_eddy
+    (rev_b0_count == 0 && rev_dwi_count == 0 && params.run_eddy) || (!params.run_topup && params.run_eddy)
 
     script:
     """
@@ -410,10 +475,10 @@ process Denoise_DWI {
     label 'big_mem'
 
     input:
-    set sid, file(dwi) from dwi_for_denoise
+    set sid, file(dwi), val(rev) from dwi_for_denoise
 
     output:
-    set sid, "${sid}__dwi_denoised.nii.gz" into\
+    set sid, "${sid}_${rev}dwi_denoised.nii.gz", val(rev) into\
         dwi_denoised_for_mix
 
     when:
@@ -427,23 +492,23 @@ process Denoise_DWI {
     export OMP_NUM_THREADS=1
     export OPENBLAS_NUM_THREADS=1
     dwidenoise $dwi dwi_denoised.nii.gz -extent $params.extent -nthreads $task.cpus
-    fslmaths dwi_denoised.nii.gz -thr 0 ${sid}__dwi_denoised.nii.gz
+    fslmaths dwi_denoised.nii.gz -thr 0 ${sid}_${rev}dwi_denoised.nii.gz
     """
 }
 
 dwi_for_test_denoise
     .map{it -> if(!params.run_dwi_denoising){it}}
     .mix(dwi_denoised_for_mix)
-    .into{dwi_for_gibbs; dwi_for_eddy; dwi_for_topup; dwi_for_eddy_topup; dwi_for_test_gibbs}
+    .into{dwi_for_gibbs; dwi_for_test_gibbs}
 
 process Gibbs_correction {
     cpus params.processes_denoise_dwi
 
     input:
-    set sid, file(dwi) from dwi_for_gibbs
+    set sid, file(dwi), val(rev) from dwi_for_gibbs
 
     output:
-    set sid, "${sid}__dwi_gibbs_corrected.nii.gz" into\
+    set sid, "${sid}_${rev}dwi_gibbs_corrected.nii.gz", val(rev) into\
         dwi_gibbs_for_mix
 
     when:
@@ -454,7 +519,7 @@ process Gibbs_correction {
     export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
     export OMP_NUM_THREADS=1
     export OPENBLAS_NUM_THREADS=1
-    mrdegibbs $dwi ${sid}__dwi_gibbs_corrected.nii.gz -nthreads $task.cpus
+    mrdegibbs $dwi ${sid}_${rev}dwi_gibbs_corrected.nii.gz -nthreads $task.cpus
     """
 }
 
@@ -463,108 +528,163 @@ dwi_for_test_gibbs
     .mix(dwi_gibbs_for_mix)
     .into{dwi_for_eddy; dwi_for_topup; dwi_for_eddy_topup; dwi_for_test_eddy_topup}
 
-dwi_for_topup
-    .join(gradients_for_topup)
-    .join(rev_b0)
-    .join(readout_encoding_for_topup)
-    .set{dwi_gradients_rev_b0_for_topup}
+sid_rev_b0_included
+    .mix(sid_rev_dwi_included_for_topup)
+    .combine(dwi_for_topup, by: 0)
+    .map{ [it[0] + it[2]] + it }
+    .join(gradients_for_prepare_topup.map{ [it[0] + it[3], it[1], it[2]] })
+    .map{ it[1..-1] }
+    .set{dwi_gradients_rev_b0_for_prepare_topup}
+
+process Prepare_for_Topup {
+  cpus 2
+
+  input:
+    set sid, file(dwi), val(rev), file(bval), file(bvec)\
+      from dwi_gradients_rev_b0_for_prepare_topup
+
+  output:
+    set sid, "${sid}_${rev}b0_mean.nii.gz", val(rev) into simple_b0_for_topup
+
+  when:
+    params.run_topup && params.run_eddy
+
+  script:
+  """
+    scil_extract_b0.py $dwi $bval $bvec ${sid}_${rev}b0_mean.nii.gz --mean\
+        --b0_thr $params.b0_thr_extract_b0 --force_b0_threshold
+  """
+}
+
+simple_b0_for_topup
+  .branch{
+    forward_b0: it[2] == "_"
+        return it[0..1]
+    reverse_b0: it[2] == "_rev_"
+        return it[0..1]
+  }
+  .set{branch_b0_for_topup}
+
+branch_b0_for_topup.reverse_b0
+  .mix(rev_b0_for_topup)
+  .join(branch_b0_for_topup.forward_b0)
+  .mix(complex_rev_b0_for_topup)
+  .join(readout_encoding_for_topup)
+  .set{rev_b0_with_readout_encoding_for_topup}
 
 process Topup {
-    cpus 2
+    cpus 4
 
     input:
-    set sid, file(dwi), file(bval), file(bvec), file(rev_b0), readout, encoding\
-        from dwi_gradients_rev_b0_for_topup
+      set sid, file(rev_b0), file(b0),  readout, encoding\
+        from rev_b0_with_readout_encoding_for_topup
 
     output:
-    set sid, "${sid}__corrected_b0s.nii.gz", "${params.prefix_topup}_fieldcoef.nii.gz",
-    "${params.prefix_topup}_movpar.txt" into topup_files_for_eddy_topup
-    file "${sid}__rev_b0_warped.nii.gz"
-    file "${sid}__b0_mean.nii.gz"
+      set sid, "${sid}__corrected_b0s.nii.gz", "${params.prefix_topup}_fieldcoef.nii.gz",
+      "${params.prefix_topup}_movpar.txt" into topup_files_for_eddy_topup
+      file "${sid}__rev_b0_warped.nii.gz"
 
     when:
-    params.run_topup && params.run_eddy
+      params.run_topup && params.run_eddy
 
     script:
     """
-    export OMP_NUM_THREADS=$task.cpus
-    export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
-    export OPENBLAS_NUM_THREADS=1
-    export ANTS_RANDOM_SEED=1234
-    scil_extract_b0.py $dwi $bval $bvec ${sid}__b0_mean.nii.gz --mean\
-        --b0_thr $params.b0_thr_extract_b0 --force_b0_threshold
-    scil_image_math.py mean $rev_b0 $rev_b0 --data_type float32 -f
-    antsRegistrationSyNQuick.sh -d 3 -f ${sid}__b0_mean.nii.gz -m $rev_b0 -o output -t r -e 1
-    mv outputWarped.nii.gz ${sid}__rev_b0_warped.nii.gz
-    scil_prepare_topup_command.py ${sid}__b0_mean.nii.gz ${sid}__rev_b0_warped.nii.gz\
-        --config $params.config_topup\
-        --encoding_direction $encoding\
-        --readout $readout --out_prefix $params.prefix_topup\
-        --out_script
-    sh topup.sh
-    cp corrected_b0s.nii.gz ${sid}__corrected_b0s.nii.gz
+      export OMP_NUM_THREADS=$task.cpus
+      export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
+      export OPENBLAS_NUM_THREADS=1
+      export ANTS_RANDOM_SEED=1234
+      antsRegistrationSyNQuick.sh -d 3 -f $b0 -m $rev_b0 -o output -t r -e 1
+      mv outputWarped.nii.gz ${sid}__rev_b0_warped.nii.gz
+      scil_prepare_topup_command.py $b0 ${sid}__rev_b0_warped.nii.gz\
+          --config $params.config_topup\
+          --encoding_direction $encoding\
+          --readout $readout --out_prefix $params.prefix_topup\
+          --out_script
+      sh topup.sh
+      cp corrected_b0s.nii.gz ${sid}__corrected_b0s.nii.gz
     """
 }
 
-dwi_for_eddy
-    .join(gradients_for_eddy)
-    .join(b0_mask_for_eddy)
-    .join(readout_encoding_for_eddy)
-    .set{dwi_gradients_mask_topup_files_for_eddy}
+dwi_for_eddy_topup.into{complex_dwi_for_eddy_topup;simple_dwi_for_eddy_topup}
+
+// Extract subjects with reverse DWI for Prepare_dwi_for_eddy
+complex_dwi_for_eddy_topup
+    .map{ [it[0] + it[2]] + it }
+    .join(gradients_for_prepare_dwi_for_eddy.map{ [it[0] + it[3], it[1], it[2]] })
+    .map{ it[1..-1] }
+    .set{dwi_gradient_for_prepare_dwi_for_eddy}
+
+sid_rev_dwi_included
+    .combine(dwi_gradient_for_prepare_dwi_for_eddy, by: 0)
+    .branch{
+        forward_dwi: it[2] == "_"
+            return it[0..1] + it[3..-1]
+        reverse_dwi: it[2] == "_rev_"
+            return it[0..1] + it[3..-1]
+    }
+    .set{branch_dwi_gradient_for_prepare_dwi_for_eddy}
 
-process Eddy {
-    cpus params.processes_eddy
+branch_dwi_gradient_for_prepare_dwi_for_eddy.forward_dwi
+    .join(branch_dwi_gradient_for_prepare_dwi_for_eddy.reverse_dwi)
+    .set{dwi_rev_gradient_for_prepare_dwi_for_eddy}
 
-    input:
-    set sid, file(dwi), file(bval), file(bvec), file(mask), readout, encoding\
-        from dwi_gradients_mask_topup_files_for_eddy
-    val(rev_b0_count) from rev_b0_counter
+// Extract subjects with reverse b0 images for Eddy
+expl1 = Channel.value(0)
 
-    output:
-    set sid, "${sid}__dwi_corrected.nii.gz" into\
-        dwi_from_eddy
-    set sid, "${sid}__bval_eddy", "${sid}__dwi_eddy_corrected.bvec" into\
-        gradients_from_eddy
+gradients_for_eddy_topup
+    .filter{ it[3] == "_" }
+    .map{ it[0..2] }
+    .set{simple_gradients_for_eddy_topup}
 
-    when:
-    (rev_b0_count == 0 && params.run_eddy) || (!params.run_topup && params.run_eddy)
+sid_rev_b0_included_for_eddy_topup
+    .combine(simple_dwi_for_eddy_topup, by: 0)
+    .filter{ it[2] == "_" }
+    .map{ it[0..1] }
+    .join(simple_gradients_for_eddy_topup)
+    .merge(expl1)
+    .set{simple_dwi_gradients_for_eddy_topup}
 
-    // Corrected DWI is clipped to 0 since Eddy can introduce negative values.
-    script:
-        slice_drop_flag=""
-        if (params.use_slice_drop_correction) {
-            slice_drop_flag="--slice_drop_correction"
-        }
-        """
-        export OMP_NUM_THREADS=$task.cpus
-        export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=$task.cpus
-        export OPENBLAS_NUM_THREADS=1
-        scil_prepare_eddy_command.py $dwi $bval $bvec $mask\
-            --eddy_cmd $params.eddy_cmd --b0_thr $params.b0_thr_extract_b0\
-            --encoding_direction $encoding\
-            --readout $readout --out_script --fix_seed\
-            $slice_drop_flag
-        sh eddy.sh
-        fslmaths dwi_eddy_corrected.nii.gz -thr 0 ${sid}__dwi_corrected.nii.gz
-        mv dwi_eddy_corrected.eddy_rotated_bvecs ${sid}__dwi_eddy_corrected.bvec
-        mv $bval ${sid}__bval_eddy
-        """
+// DOES NOT WORK FOR EVERYTHING BECAUSE AP PA which one is rev, LR RL same question
+
+process Prepare_dwi_for_eddy {
+  cpus 2
+
+  input:
+    set sid, file(dwi), file(bval), file(bvec), file(rev_dwi), file(rev_bval), \
+        file(rev_bvec) from dwi_rev_gradient_for_prepare_dwi_for_eddy
+
+  output:
+    set sid, "${sid}__concatenated_dwi.nii.gz", "${sid}__concatenated_dwi.bval", "${sid}__concatenated_dwi.bvec", env(rev_number_dir) into concatenated_dwi_for_eddy
+
+  when:
+    params.run_topup && params.run_eddy
+
+  script:
+  """
+    scil_concatenate_dwi.py ${sid}__concatenated_dwi.nii.gz ${sid}__concatenated_dwi.bval ${sid}__concatenated_dwi.bvec -f\
+      --in_dwis ${dwi} ${rev_dwi} --in_bvals ${bval} ${rev_bval}\
+      --in_bvecs ${bvec} ${rev_bvec}
+
+    rev_number_dir=\$(scil_print_header.py ${rev_dwi} --key dim | sed "s/  / /g" | sed "s/  / /g" | rev | cut -d' ' -f4-4 | rev)
+  """
 }
 
-dwi_for_eddy_topup
-    .join(gradients_for_eddy_topup)
+concatenated_dwi_for_eddy
+    .mix(simple_dwi_gradients_for_eddy_topup)
     .join(topup_files_for_eddy_topup)
     .join(readout_encoding_for_eddy_topup)
     .set{dwi_gradients_mask_topup_files_for_eddy_topup}
 
 process Eddy_Topup {
     cpus params.processes_eddy
+    memory { 5.GB * task.attempt }
 
     input:
-    set sid, file(dwi), file(bval), file(bvec), file(b0s_corrected),
+    set sid, file(dwi), file(bval), file(bvec), val(number_rev_dwi), file(b0s_corrected),
         file(field), file(movpar), readout, encoding\
         from dwi_gradients_mask_topup_files_for_eddy_topup
     val(rev_b0_count) from rev_b0_counter
+    val(rev_dwi_count) from rev_dwi_counter
 
     output:
     set sid, "${sid}__dwi_corrected.nii.gz" into\
@@ -574,7 +694,7 @@ process Eddy_Topup {
     file "${sid}__b0_bet_mask.nii.gz"
 
     when:
-    rev_b0_count > 0 && params.run_topup && params.run_eddy
+    (rev_b0_count > 0 || rev_dwi_count > 0) && params.run_topup && params.run_eddy
 
     // Corrected DWI is clipped to ensure there are no negative values
     // introduced by Eddy.
@@ -594,6 +714,64 @@ process Eddy_Topup {
             --b0_thr $params.b0_thr_extract_b0\
             --encoding_direction $encoding\
             --readout $readout --out_script --fix_seed\
+            --n_reverse ${number_rev_dwi}\
+            --lsr_resampling\
+            $slice_drop_flag
+	echo "--very_verbose" >> eddy.sh
+	sh eddy.sh
+        fslmaths dwi_eddy_corrected.nii.gz -thr 0 ${sid}__dwi_corrected.nii.gz
+
+	if [[ $number_rev_dwi -eq 0 ]]
+	then
+	   mv dwi_eddy_corrected.eddy_rotated_bvecs ${sid}__dwi_eddy_corrected.bvec
+          mv $bval ${sid}__bval_eddy
+	else
+	   scil_validate_and_correct_eddy_gradients.py dwi_eddy_corrected.eddy_rotated_bvecs $bval ${number_rev_dwi} ${sid}__dwi_eddy_corrected.bvec ${sid}__bval_eddy
+	fi
+	"""
+}
+
+dwi_for_eddy
+    .map{ [it[0] + it[2]] + it }
+    .join(gradients_for_eddy.map{ [it[0] + it[3], it[1], it[2]] })
+    .filter{ it[3] == "_" }
+    .map{ it[1..2] + it[4..-1] }
+    .join(b0_mask_for_eddy)
+    .join(readout_encoding_for_eddy)
+    .set{dwi_gradients_mask_topup_files_for_eddy}
+
+process Eddy {
+    cpus params.processes_eddy
+
+    input:
+    set sid, file(dwi), file(bval), file(bvec), file(mask), readout, encoding\
+        from dwi_gradients_mask_topup_files_for_eddy
+    val(rev_b0_count) from rev_b0_counter
+    val(rev_dwi_count) from rev_dwi_counter
+
+    output:
+    set sid, "${sid}__dwi_corrected.nii.gz" into\
+        dwi_from_eddy
+    set sid, "${sid}__bval_eddy", "${sid}__dwi_eddy_corrected.bvec" into\
+        gradients_from_eddy
+
+    when:
+    (rev_b0_count == 0 && rev_dwi_count == 0 && params.run_eddy) || (!params.run_topup && params.run_eddy)
+
+    // Corrected DWI is clipped to 0 since Eddy can introduce negative values.
+    script:
+        slice_drop_flag=""
+        if (params.use_slice_drop_correction) {
+            slice_drop_flag="--slice_drop_correction"
+        }
+        """
+        export OMP_NUM_THREADS=$task.cpus
+        export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=$task.cpus
+        export OPENBLAS_NUM_THREADS=1
+        scil_prepare_eddy_command.py $dwi $bval $bvec $mask\
+            --eddy_cmd $params.eddy_cmd --b0_thr $params.b0_thr_extract_b0\
+            --encoding_direction $encoding\
+            --readout $readout --out_script --fix_seed\
             $slice_drop_flag
         sh eddy.sh
         fslmaths dwi_eddy_corrected.nii.gz -thr 0 ${sid}__dwi_corrected.nii.gz
@@ -602,8 +780,18 @@ process Eddy_Topup {
         """
 }
 
-dwi_for_test_eddy_topup.map{it -> if(!params.run_eddy){it}}.set{dwi_for_skip_eddy_topup}
-gradients_for_test_eddy_topup.map{it -> if(!params.run_eddy){it}}.set{gradients_for_skip_eddy_topup}
+
+dwi_for_test_eddy_topup
+    .map{it -> if(!params.run_eddy){it}}
+    .filter{ it[2] == "_" }
+    .map{ it[0..1] }
+    .set{dwi_for_skip_eddy_topup}
+
+gradients_for_test_eddy_topup
+    .map{it -> if(!params.run_eddy){it}}
+    .filter{ it[3] == "_" }
+    .map{ it[0..2] }
+    .set{gradients_for_skip_eddy_topup}
 
 dwi_from_eddy
     .mix(dwi_from_eddy_topup)
diff --git a/nextflow.config b/nextflow.config
old mode 100755
new mode 100644
index edcf3a4..632d0d9
--- a/nextflow.config
+++ b/nextflow.config
@@ -193,6 +193,7 @@ profiles {
 
     use_cuda {
         singularity.runOptions='--nv'
+        docker.runOptions='--gpus all'
         params.eddy_cmd="eddy_cuda"
     }
 
@@ -262,5 +263,4 @@ profiles {
         params.pft_nbr_seeds=30
         params.pft_seeding_mask_type="interface"
     }
-
 }