DWI and SC updates: Complex data denoising

functions/02_proc-dwi.sh

@@ -23,16 +23,17 @@ nocleanup=$5
 threads=$6
 tmpDir=$7
 dwi_main=$8
-dwi_rpe=$9
-dwi_processed=${10}
-rpe_all=${11}
-regAffine=${12}
-dwi_str=${13}
-b0thr=${14}
-bvalscale=${15}
-synth_reg=${16}
-dwi_upsample=${17}
-PROC=${18}
+dwi_phase=$9
+dwi_rpe=${10}
+dwi_processed=${11}
+rpe_all=${12}
+regAffine=${13}
+dwi_str=${14}
+b0thr=${15}
+bvalscale=${16}
+synth_reg=${17}
+dwi_upsample=${18}
+PROC=${19}
 here=$(pwd)
 
 #------------------------------------------------------------------------------#
@@ -51,6 +52,7 @@ bids_variables "$BIDS" "$id" "$out" "$SES"
 Info "Inputs of proc_dwi"
 Note "tmpDir        :" "$tmpDir"
 Note "dwi_main      :" "$dwi_main"
+Note "dwi_phase     :" "$dwi_phase"
 Note "dwi_rpe       :" "$dwi_rpe"
 Note "rpe_all       :" "$rpe_all"
 Note "dwi_acq       :" "$dwi_str"
@@ -80,6 +82,13 @@ if [[ "$dwi_main" != "DEFAULT" ]]; then
     IFS=',' read -ra bids_dwis <<< "$dwi_main"
     bids_dwis=("${bids_dwis[@]}")
 fi
+# Manage manual inputs: DWI phase image(s)
+if [[ "$dwi_phase" != "DEFAULT" && "$dwi_phase" != "FALSE" ]]; then
+    IFS=',' read -ra bids_dwis_phase <<< "$dwi_phase"
+    bids_dwis_phase=("${bids_dwis_phase[@]}")
+elif [[ "$dwi_phase" == "FALSE" ]]; then
+    bids_dwis_phase=()
+fi
 # Manage manual inputs: DWI reverse phase encoding
 if [[ "$dwi_rpe" != "DEFAULT" ]]; then
     IFS=',' read -ra dwi_reverse <<< "$dwi_rpe"
@@ -165,50 +174,55 @@ b0_refacq=$(echo "${bids_dwis[0]##*/}" | awk -F ".nii" '{print $1}'); b0_refacq=
 if [[ "$dwi_processed" == "FALSE" ]] && [[ ! -f "$dwi_corr" ]]; then
     if [ ! -f "$dwi_res" ] || [ ! -f "$dwi_dns" ]; then ((N++))
     Info "DWI denoise and concatenation"
-    # Concatenate shells -if only one shell then just convert to mif and rename.
-          for dwi in "${bids_dwis[@]}"; do
+    # Denoise and remove gibbs ringing artifact.
+          for idx in "${!bids_dwis[@]}"; do
+                dwi="${bids_dwis[$idx]}"
                 dwi_nom=$(echo "${dwi##*/}" | awk -F ".nii" '{print $1}')
-                bids_dwi_str=$(echo "$dwi" | awk -F . '{print $1}')
-                Do_cmd mrconvert "$dwi" -json_import "${bids_dwi_str}.json" -fslgrad "${bids_dwi_str}.bvec" "${bids_dwi_str}.bval" "${tmp}/${dwi_nom}.mif" "${bvalstr}"
-                Do_cmd dwiextract "${tmp}/${dwi_nom}.mif" "${tmp}/${dwi_nom}_b0.mif" -bzero
-                Do_cmd mrmath "${tmp}/${dwi_nom}_b0.mif" mean "${tmp}/${dwi_nom}_b0.nii.gz" -axis 3 -nthreads "$threads"
+                bids_dwi_str=$(echo "$dwi" | awk -F "_dwi." '{print $1}')
+                Do_cmd mrconvert "$dwi" -json_import "${bids_dwi_str}_dwi.json" -fslgrad "${bids_dwi_str}_dwi.bvec" "${bids_dwi_str}_dwi.bval" "${tmp}/${dwi_nom}.mif" "${bvalstr}"
+                # lookup if dwi has phase images
+                if [[ -n "${bids_dwis_phase[$idx]:-}" ]]; then
+                    phase_dwi="${bids_dwis_phase[$idx]}"
+                    phase_nom=$(basename "$phase_dwi" | awk -F ".nii" '{print $1}')
+                    phase_str=$(echo "$phase_dwi" | awk -F "_dwi." '{print $1}')
+                    Do_cmd mrconvert "$phase_dwi" -json_import "${phase_str}_dwi.json" -fslgrad "${phase_str}_dwi.bvec" "${phase_str}_dwi.bval" "${tmp}/${phase_nom}.mif" "${bvalstr}"
+                    Do_cmd designer -denoise -shrinkage frob -algorithm jespersen -phase "${tmp}/${phase_nom}.mif" -degibbs -scratch ${tmp}/designer_proc "${tmp}/${dwi_nom}.mif" "${tmp}/${dwi_nom}_dn_tmp.mif"
+                    Do_cmd mrconvert "${tmp}/${dwi_nom}_dn_tmp.mif" "${tmp}/${dwi_nom}_dn.mif" -json_import "${bids_dwi_str}_dwi.json"
+                else
+                    Do_cmd designer -denoise -shrinkage frob -adaptive_patch -rician -degibbs -scratch ${tmp}/designer_proc "${tmp}/${dwi_nom}.mif" "${tmp}/${dwi_nom}_dn_tmp.mif"
+                    Do_cmd mrconvert "${tmp}/${dwi_nom}_dn_tmp.mif" "${tmp}/${dwi_nom}_dn.mif"  -json_import "${bids_dwi_str}_dwi.json"
+                fi
+                Do_cmd dwiextract "${tmp}/${dwi_nom}_dn.mif" "${tmp}/${dwi_nom}_dn_b0.mif" -bzero
+                Do_cmd mrmath "${tmp}/${dwi_nom}_dn_b0.mif" mean "${tmp}/${dwi_nom}_dn_b0.nii.gz" -axis 3 -nthreads "$threads"
           done
 
-          # Rigid registration between shells
+          # Rigid registration between input images -- if multiple files
           n=$((${#bids_dwis[*]} - 1))
           if [[ ${#bids_dwis[*]} -gt 1 ]]; then
-            b0_ref=${tmp}/$(echo "${bids_dwis[0]##*/}" | awk -F ".nii" '{print $1}')_b0.nii.gz
+            b0_ref=${tmp}/$(echo "${bids_dwis[0]##*/}" | awk -F ".nii" '{print $1}')_dn_b0.nii.gz
             for ((i=1; i<=n; i++)); do
                 dwi_nom=$(echo "${bids_dwis[i]##*/}" | awk -F ".nii" '{print $1}')
                 bids_dwi_str=$(echo "${bids_dwis[i]}" | awk -F . '{print $1}')
-                b0_nom="${tmp}/$(echo "${bids_dwis[i]##*/}" | awk -F ".nii" '{print $1}')_b0.nii.gz"
+                b0_nom="${tmp}/$(echo "${bids_dwis[i]##*/}" | awk -F ".nii" '{print $1}')_dn_b0.nii.gz"
                 b0_acq=$(echo "${dwi_nom/${idBIDS}_/}"); b0_acq=$(echo "${b0_acq/_dwi/}")
                 b0mat_str="${tmp}/${idBIDS}_from-${b0_acq}_to-${b0_refacq}${dwi_str_}_mode-image_desc-rigid_"
                 b0mat="${b0mat_str}0GenericAffine.mat"
 
                 Info "Registering ${b0_acq} to ${b0_refacq}"
                 Do_cmd antsRegistrationSyN.sh -d 3 -m "$b0_nom" -f "$b0_ref" -o "$b0mat_str" -t r -n "$threads" -p d
-                mrconvert "${tmp}/${dwi_nom}.mif" "${tmp}/${dwi_nom}.nii.gz"
-                Do_cmd antsApplyTransforms -d 3 -e 3 -i "${tmp}/${dwi_nom}.nii.gz" -r "$b0_ref" -t "$b0mat" -o "${tmp}/${dwi_nom}_in-${b0_refacq}.nii.gz" -v -u int
-                Do_cmd mrconvert "${tmp}/${dwi_nom}_in-${b0_refacq}.nii.gz" -json_import "${bids_dwi_str}.json" -fslgrad "${bids_dwi_str}.bvec" "${bids_dwi_str}.bval" "${tmp}/${dwi_nom}__Ralign.mif" -force -quiet  "${bvalstr}"
+                Do_cmd ConvertTransformFile 3 "${b0mat}" "${tmp}/ants.txt"
+                Do_cmd transformconvert "${tmp}/ants.txt" itk_import "${tmp}/mrtrix_affine.txt"
+                Do_cmd mrtransform "${tmp}/${dwi_nom}_dn.mif" "${tmp}/${dwi_nom}_dn__Ralign.mif" -linear "${tmp}/mrtrix_affine.txt"
             done
           fi
 
           Info "Concatenatenating shells"
           dwi_0=$(echo "${bids_dwis[0]##*/}" | awk -F ".nii" '{print $1}')
           if [ "${#bids_dwis[@]}" -eq 1 ]; then
-            cp "${tmp}/${dwi_0}.mif" "$dwi_cat"
+            cp "${tmp}/${dwi_0}_dn.mif" "$dwi_dns"
           else
-            Do_cmd mrcat "${tmp}/${dwi_0}.mif" "${tmp}/*_Ralign.mif" "$dwi_cat" -nthreads "$threads"
+            Do_cmd mrcat "${tmp}/${dwi_0}_dn.mif" "${tmp}/*_dn__Ralign.mif" "$dwi_dns" -nthreads "$threads"
           fi
-
-          # Denoise DWI and calculate residuals
-          Info "DWI MP-PCA denoising and Gibbs ringing correction"
-          dwi_dns_tmp="${tmp}/MP-PCA_dwi.mif"
-          Do_cmd dwidenoise "$dwi_cat" "$dwi_dns_tmp" -nthreads "$threads"
-          mrcalc "$dwi_cat" "$dwi_dns_tmp" -subtract - -nthreads "$threads" | mrmath - mean "$dwi_resPCA" -axis 3
-          Do_cmd mrdegibbs "$dwi_dns_tmp" "$dwi_dns" -nthreads "$threads"
-          mrcalc "$dwi_dns_tmp" "$dwi_dns" -subtract - -nthreads "$threads" | mrmath - mean "$dwi_resGibss" -axis 3
           ((Nsteps++))
     else
           Info "Subject ${id} has DWI in mif, denoised and concatenaded"; ((Nsteps++)); ((N++))
@@ -229,7 +243,13 @@ if [[ "$dwi_processed" == "FALSE" ]] && [[ ! -f "$dwi_corr" ]]; then
                 Info "DWI reverse phase encoding processing: $dwi_nom"
                 if [[ -f "${bids_dwi_str}.bvec" ]] && [[ -f "${bids_dwi_str}.bval" ]]; then
                     Do_cmd mrconvert "$dwi" -json_import "${bids_dwi_str}.json" -fslgrad "${bids_dwi_str}.bvec" "${bids_dwi_str}.bval" "${tmp}/${dwi_nom}.mif" "${bvalstr}"
-                    Do_cmd dwiextract "${tmp}/${dwi_nom}.mif" "${tmp}/${dwi_nom}_b0.mif" -bzero
+                    rpe_dim=$(mrinfo "${tmp}/${dwi_nom}.mif" -ndim)
+                    if [[ "$rpe_dim" -eq 3 ]]; then
+                      Warning "Only on rpe volume was found the script will assumme that the volume is b0!!!!"
+                      Do_cmd cp "${tmp}/${dwi_nom}.mif" "${tmp}/${dwi_nom}_b0.mif"
+                    elif [[ "$rpe_dim" -gt 3 ]]; then
+                      Do_cmd dwiextract "${tmp}/${dwi_nom}.mif" "${tmp}/${dwi_nom}_b0.mif" -bzero
+                    fi
                 else
                     Warning "No bval or bvecs were found the script will assumme that all the volumes are b0s!!!!"
                     Do_cmd mrconvert "$dwi" -json_import "${bids_dwi_str}.json" "${tmp}/${dwi_nom}.mif" "${bvalstr}"
@@ -370,7 +390,7 @@ if [[ ! -f "$dwi_corr" ]]; then ((N++))
       # DWIs all acquired with a single fixed phase encoding; but additionally a
       # pair of b=0 images with reversed phase encoding to estimate the inhomogeneity field:
       echo -e "COMMAND --> dwifslpreproc $dwi_4proc $dwi_corr $opt -pe_dir $pe_dir -readout_time $ReadoutTime -eddy_options \" --data_is_shelled --slm=linear --repol\" -nthreads $threads -nocleanup -scratch $tmp -force"
-      dwifslpreproc "$dwi_4proc" "$dwi_corr" $opt -pe_dir "$pe_dir" -readout_time "$ReadoutTime" -eddy_options " --data_is_shelled --slm=linear --repol" -nthreads "$threads" -nocleanup -scratch "$tmp" -force
+      dwifslpreproc "$dwi_4proc" "$dwi_corr" $opt -pe_dir "$pe_dir" -readout_time "$ReadoutTime" -eddy_options " --data_is_shelled --slm=linear --repol" -eddyqc_text "${tmp}/eddy_qc" -nthreads "$threads" -nocleanup -scratch "$tmp" -force
       # Step QC
       if [[ ! -f "$dwi_corr" ]]; then Error "dwifslpreproc failed, check the logs"; exit;
       else

functions/utilities.sh

@@ -95,7 +95,8 @@ export idBIDS="${subject}${ses}"
 
   # BIDS Files
   bids_T1ws=($(ls "$subject_bids"/anat/*T1w.nii* 2>/dev/null))
-  bids_dwis=($(ls "${subject_bids}/dwi/${subject}${ses}"*_dir-AP_*dwi.nii* 2>/dev/null))
+  bids_dwis=($(ls "${subject_bids}/dwi/${subject}${ses}"*_dir-AP_dwi.nii* 2>/dev/null))
+  bids_dwis_phase=($(ls "${subject_bids}/dwi/${subject}${ses}"*_dir-AP_part-phase_dwi.nii* 2>/dev/null))
   bids_T1map=$(ls "$subject_bids"/anat/*mp2rage*T1map.nii* 2>/dev/null)
   bids_inv1=$(ls "$subject_bids"/anat/*inv1*T1map.nii* 2>/dev/null)
   bids_inv2=$(ls "$subject_bids"/anat/*inv2*T1map.nii* 2>/dev/null)
@@ -185,7 +186,12 @@ bids_print.variables-dwi() {
   # This functions prints BIDS variables names and files if found
   Info "Variables for DWI processing"
   Note "proc_dwi dir    :" "$proc_dwi"
-  Note "bids_dwis       :" "N-${#bids_dwis[@]}, $bids_dwis"
+  for i in "${!bids_dwis[@]}"; do
+    file.exist "bids_dwis $[$i+1]/${#bids_dwis[@]}       :" ${bids_dwis[i]}
+  done
+  for i in "${!bids_dwis_phase[@]}"; do
+    file.exist "bids_dwis_phase $[$i+1]/${#bids_dwis[@]} :" ${bids_dwis_phase[i]}
+  done
   Note "dwi_reverse     :" "N-${#dwi_reverse[@]}, $dwi_reverse"
 
   Note "T1 nativepro    :" "$(find "$T1nativepro" 2>/dev/null)"
@@ -261,6 +267,7 @@ bids_variables_unset() {
   unset icafixTraining
   unset bids_T1ws
   unset bids_dwis
+  unset bids_dwis_phase
   unset bids_T1map
   unset bids_inv1
   unset dwi_reverse

micapipe

@@ -82,6 +82,9 @@ echo -e "
 \t\033[38;5;197m-proc_dwi\033[0m
 \t   \033[38;5;120m-dwi_main\033[0m            : Full path to DWI with N number of directions and b0. If used, it will overwrite the defaults.
 \t\t\t            Default = <bids>/<subject>/dwi/*_dir-AP_dwi.nii*
+\t   \033[38;5;120m-dwi_phase\033[0m           : Full path to DWI phase part with N number of directions. If used, it will overwrite the defaults.
+\t\t\t            Default = <bids>/<subject>/dwi/*_dir-AP_part-phase_dwi.nii*
+\t\t\t            Set to 'FALSE' if no DWI phase images will be provided
 \t   \033[38;5;120m-dwi_rpe\033[0m             : DWI b0 image(s) with phase reversal encoding. If used, it will overwrite the defaults.
 \t\t\t            Default = <bids>/<subject>/dwi/<sub>_dir-PA_dwi.nii*
 \t\t\t            Set to 'FALSE' if no DWI with reverse fase encoding will be provided
@@ -455,6 +458,10 @@ do
     dwi_main="$2"
     shift;shift;
   ;;
+  -dwi_phase)
+    dwi_phase="$2"
+    shift;shift;
+  ;;
   -dwi_rpe)
     dwi_rpe="$2"
     shift;shift;
@@ -608,6 +615,16 @@ else
     for i in "${!dwis_main[@]}"; do dwis_main[i]=$(realpath ${dwis_main[$i]}); done # Full path
     dwi_main=$(IFS=','; echo "${dwis_main[*]}")
 fi
+# User provided DWI phase images
+if [ -z ${dwi_phase} ]; then
+    dwi_phase=DEFAULT
+elif [ ${dwi_phase} == "FALSE" ]; then
+    dwi_phase=${dwi_phase}
+else
+    IFS=',' read -ra dwis_phase <<< $dwi_phase
+    for i in "${!dwis_phase[@]}"; do dwis_phase[i]=$(realpath ${dwis_phase[$i]}); done # Full path
+    dwi_phase=$(IFS=','; echo "${dwis_phase[*]}")
+fi
 # User provided DWI reverse phase encoded
 if [ -z ${dwi_rpe} ]; then
     dwi_rpe=DEFAULT
@@ -832,7 +849,7 @@ fi
 if [ "$procDWI" = "TRUE" ]; then
     if [[ ${dwi_str} != "DEFAULT" ]]; then log_str=_${dwi_str}; else log_str=""; fi
     log_file_str=$dir_logs/proc_dwi_$(date +'%Y-%m-%d_%H.%M.%S')${log_str}
-    COMMAND="${scriptDir}/02_proc-dwi.sh $BIDS $id $out $SES $nocleanup $threads $tmpDir $dwi_main $dwi_rpe $dwi_processed $rpe_all $regAffine $dwi_str $b0thr $bvalscale $synth_reg $dwi_upsample"
+    COMMAND="${scriptDir}/02_proc-dwi.sh $BIDS $id $out $SES $nocleanup $threads $tmpDir $dwi_main $dwi_phase $dwi_rpe $dwi_processed $rpe_all $regAffine $dwi_str $b0thr $bvalscale $synth_reg $dwi_upsample"
     # mica.q - Diffusion processing
     if [[ $micaq == "TRUE" ]]; then
         Info "MICA qsub - Diffusion processing"