ENH: Generate the eddy QC report

DiffusionPreprocessing/DiffPreprocPipeline_PostEddy.sh

@@ -353,6 +353,10 @@ main()
 		cp -p ${filename} ${to_location}
 	done
 
+  mkdir -p ${outdirT1w/qc
+  cp ${outdir}/qc/* ${outdirT1w}/qc
+  immv ${outdirT1w}/cnr_maps ${outdirT1w}/qc/cnr_maps
+
 	log_Msg "Completed!"
 	exit 0
 }

DiffusionPreprocessing/scripts/DiffusionToStructural.sh

@@ -106,24 +106,34 @@ if [ ${GdcorrectionFlag} -eq 1 ]; then
     # Dilation outside of the field of view to minimise the effect of the hard field of view edge on the interpolation
     ${CARET7DIR}/wb_command -volume-dilate $DataDirectory/warped/data_warped.nii.gz $DilateDistance NEAREST $DataDirectory/warped/data_warped_dilated.nii.gz
     ${FSLDIR}/bin/applywarp --rel -i "$DataDirectory"/warped/data_warped_dilated -r "$T1wRestoreImage"_${DiffRes} -w "$WorkingDirectory"/grad_unwarp_diff2str --interp=spline -o "$T1wOutputDirectory"/data
+    ${FSLDIR}/bin/imrm $DataDirectory/warped/data_warped_dilated
 
     # Transforms field of view mask to T1-weighted space
     # (Be sure to use the fov_mask derived prior to application of GDC)
     ${FSLDIR}/bin/applywarp --rel -i "$DataDirectory"/warped/fov_mask_warped -r "$T1wRestoreImage"_${DiffRes} -w "$WorkingDirectory"/grad_unwarp_diff2str --interp=trilinear -o "$T1wOutputDirectory"/fov_mask
 
+    # Transform CNR maps
+    ${CARET7DIR}/wb_command -volume-dilate $DataDirectory/warped/cnr_maps_warped.nii.gz $DilateDistance NEAREST $DataDirectory/warped/cnr_maps_dilated.nii.gz
+    ${FSLDIR}/bin/applywarp --rel -i "$DataDirectory"/warped/cnr_maps_dilated -r "$T1wRestoreImage"_${DiffRes} -w "$WorkingDirectory"/grad_unwarp_diff2str --interp=spline -o "$T1wOutputDirectory"/cnr_maps
+    ${FSLDIR}/bin/imrm $DataDirectory/warped/cnr_maps_dilated
+
     # Now register the grad_dev tensor
     ${FSLDIR}/bin/vecreg -i "$DataDirectory"/grad_dev -o "$T1wOutputDirectory"/grad_dev -r "$T1wRestoreImage"_${DiffRes} -t "$WorkingDirectory"/diff2str.mat --interp=spline
     ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/grad_dev -mas "$T1wOutputDirectory"/nodif_brain_mask_temp "$T1wOutputDirectory"/grad_dev  #Mask-out values outside the brain 
-    ${FSLDIR}/bin/imrm $DataDirectory/warped/data_warped_dilated
 else
     # Dilation outside of the field of view to minimise the effect of the hard field of view edge on the interpolation
     ${CARET7DIR}/wb_command -volume-dilate $DataDirectory/data.nii.gz $DilateDistance NEAREST $DataDirectory/data_dilated.nii.gz
     # Register diffusion data to T1w space without considering gradient nonlinearities
     ${FSLDIR}/bin/flirt -in "$DataDirectory"/data_dilated -ref "$T1wRestoreImage"_${DiffRes} -applyxfm -init "$WorkingDirectory"/diff2str.mat -interp spline -out "$T1wOutputDirectory"/data
+    ${FSLDIR}/bin/imrm $DataDirectory/data_dilated
+
+    # Transform CNR maps
+    ${CARET7DIR}/wb_command -volume-dilate $DataDirectory/cnr_maps.nii.gz $DilateDistance NEAREST $DataDirectory/cnr_maps_dilated.nii.gz
+    ${FSLDIR}/bin/flirt -in "$DataDirectory"/cnr_maps_dilated -ref "$T1wRestoreImage"_${DiffRes} -applyxfm -init "$WorkingDirectory"/diff2str.mat -interp spline -out "$T1wOutputDirectory"/cnr_maps
+    ${FSLDIR}/bin/imrm $DataDirectory/cnr_maps_dilated
 
     # Transforms field of view mask to T1-weighted space
     ${FSLDIR}/bin/flirt -in "$DataDirectory"/fov_mask -ref "$T1wRestoreImage"_${DiffRes} -applyxfm -init "$WorkingDirectory"/diff2str.mat -interp trilinear -out "$T1wOutputDirectory"/fov_mask
-    ${FSLDIR}/bin/imrm $DataDirectory/data_dilated
 fi
 
 # only include voxels fully(!) within the field of view for every volume
@@ -134,6 +144,8 @@ ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/data -mas "$T1wOutputDirectory"/nod
 ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/data -thr 0 "$T1wOutputDirectory"/data      #Remove negative intensity values (from eddy) from final data
 ${FSLDIR}/bin/imrm "$T1wOutputDirectory"/nodif_brain_mask_temp
 
+${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/cnr_maps -mas "$T1wOutputDirectory"/fov_mask "$T1wOutputDirectory"/cnr_maps
+
 # Identify any voxels that are zeros across all frames and remove those from the final nodif_brain_mask
 ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/data -Tmean "$T1wOutputDirectory"/temp
 ${FSLDIR}/bin/immv "$T1wOutputDirectory"/nodif_brain_mask.nii.gz "$T1wOutputDirectory"/nodif_brain_mask_old.nii.gz

DiffusionPreprocessing/scripts/eddy_postproc.sh

@@ -17,6 +17,20 @@ globalscriptsdir=${HCPPIPEDIR_Global}
 eddydir=${workingdir}/eddy
 datadir=${workingdir}/data
 
+echo "Generating eddy QC report in ${workingdir}/qc"
+if [ -d "{workingdir}/qc" ]; then rm -r ${workingdir}/qc; fi
+qc_command=("${FSLDIR}/bin/eddy_quad")
+qc_command+=("${eddydir}/eddy_unwarped_images")
+qc_command+=(-idx "${eddydir}/index.txt")
+qc_command+=(-par "${eddydir}/acqparams.txt")
+qc_command+=(-m "${eddydir}/nodif_brain_mask.nii.gz")
+qc_command+=(-b "${eddydir}/Pos_Neg.bvals")
+qc_command+=(-g "${eddydir}/eddy_unwarped_images.eddy_rotated_bvecs")
+qc_command+=(-o "${workingdir}/qc")
+qc_command+=(-f "${workingdir}/topup/topup_Pos_Neg_b0_field.nii.gz")
+qc_command+=(-v)
+${qc_command[@]}
+
 #Prepare for next eddy Release
 #if [ ! -e ${eddydir}/${EddyJacFlag} ]; then 
 #    echo "LSR resampling has been used. Eddy Output has already been combined."
@@ -27,7 +41,8 @@ datadir=${workingdir}/data
 if [ ${CombineDataFlag} -eq 2 ]; then
 	${FSLDIR}/bin/imcp  ${eddydir}/eddy_unwarped_images ${datadir}/data
 	cp ${eddydir}/Pos_Neg.bvals ${datadir}/bvals
-	cp ${eddydir}/Pos_Neg.bvecs ${datadir}/bvecs
+	cp ${datadir}/bvecs ${datadir}/bvecs_noRot
+	cp ${eddydir}/eddy_unwarped_images.eddy_rotated_bvecs ${datadir}/bvecs
 else
 	echo "JAC resampling has been used. Eddy Output is now combined."
 	PosVols=`wc ${eddydir}/Pos.bval | awk {'print $2'}`
@@ -87,6 +102,7 @@ else
 fi
 #fi
 
+imcp ${eddydir}/eddy_unwarped_images.eddy_cnr_maps ${datadir}/cnr_maps
 
 # Create a mask representing voxels within the field of view for all volumes prior to dilation
 # 'eddy' can return negative values in some low signal locations, so use -abs for determining the fov mask
@@ -95,42 +111,45 @@ ${FSLDIR}/bin/fslmaths ${datadir}/data -abs -Tmin -bin -fillh ${datadir}/fov_mas
 
 if [ ! $GdCoeffs = "NONE" ] ; then
     echo "Correcting for gradient nonlinearities"
-	# Note: "data_warped" is eddy-current and suspectibility distortion corrected (via 'eddy'), but prior to gradient distortion correction
-	# i.e., "data_posteddy_preGDC" would be another way to think of it
-    ${FSLDIR}/bin/immv ${datadir}/data ${datadir}/data_warped
+	  # Note: data in the warped directory is eddy-current and suspectibility distortion corrected (via 'eddy'), but prior to gradient distortion correction
+	  # i.e., "data_posteddy_preGDC" would be another way to think of it
+    mkdir -p ${datadir}/warped
+    ${FSLDIR}/bin/immv ${datadir}/data ${datadir}/warped/data_warped
+    ${FSLDIR}/bin/immv ${datadir}/fov_mask ${datadir}/warped/fov_mask_warped
+    ${FSLDIR}/bin/immv ${datadir}/cnr_maps ${datadir}/warped/cnr_maps_warped
 
     # Dilation outside of the field of view to minimise the effect of the hard field of view edge on the interpolation
-	DiffRes=`${FSLDIR}/bin/fslval ${datadir}/data_warped pixdim1`
+	  DiffRes=`${FSLDIR}/bin/fslval ${datadir}/warped/data_warped pixdim1`
     DilateDistance=`echo "$DiffRes * 4" | bc`  # Extrapolates the diffusion data up to 4 voxels outside of the FOV
-    ${CARET7DIR}/wb_command -volume-dilate ${datadir}/data_warped.nii.gz $DilateDistance NEAREST ${datadir}/data_warped_dilated.nii.gz
+    ${CARET7DIR}/wb_command -volume-dilate ${datadir}/warped/data_warped.nii.gz $DilateDistance NEAREST ${datadir}/warped/data_dilated.nii.gz
 
     # apply gradient distortion correction
-    ${globalscriptsdir}/GradientDistortionUnwarp.sh --workingdir="${datadir}" --coeffs="${GdCoeffs}" --in="${datadir}/data_warped_dilated" --out="${datadir}/data" --owarp="${datadir}/fullWarp"
+    ${globalscriptsdir}/GradientDistortionUnwarp.sh --workingdir="${datadir}" --coeffs="${GdCoeffs}" --in="${datadir}/warped/data_dilated" --out="${datadir}/data" --owarp="${datadir}/fullWarp"
+    ${FSLDIR}/bin/immv ${datadir}/fullWarp ${datadir}/warped
+    ${FSLDIR}/bin/immv ${datadir}/fullWarp_abs ${datadir}/warped
+    ${FSLDIR}/bin/imrm ${datadir}/warped/data_dilated
+
+    # Transform CNR maps
+    ${CARET7DIR}/wb_command -volume-dilate ${datadir}/warped/cnr_maps_warped.nii.gz $DilateDistance NEAREST ${datadir}/warped/cnr_maps_dilated.nii.gz
+    ${FSLDIR}/bin/applywarp --rel --interp=trilinear -i ${datadir}/warped/cnr_maps_dilated -r ${datadir}/warped/cnr_maps_dilated -w ${datadir}/warped/fullWarp -o ${datadir}/cnr_maps
+    ${FSLDIR}/bin/imrm ${datadir}/warped/cnr_maps_dilated
 
     # Transform field of view mask (using conservative trilinear interpolation with high threshold)
-    ${FSLDIR}/bin/immv ${datadir}/fov_mask ${datadir}/fov_mask_warped
-    ${FSLDIR}/bin/applywarp --rel --interp=trilinear -i ${datadir}/fov_mask_warped -r ${datadir}/fov_mask_warped -w ${datadir}/fullWarp -o ${datadir}/fov_mask
+    ${FSLDIR}/bin/applywarp --rel --interp=trilinear -i ${datadir}/warped/fov_mask_warped -r ${datadir}/warped/fov_mask_warped -w ${datadir}/warped/fullWarp -o ${datadir}/fov_mask
     ${FSLDIR}/bin/fslmaths ${datadir}/fov_mask -thr 0.999 -bin ${datadir}/fov_mask
 
     echo "Computing gradient coil tensor to correct for gradient nonlinearities"
-    ${FSLDIR}/bin/calc_grad_perc_dev --fullwarp=${datadir}/fullWarp -o ${datadir}/grad_dev
+    ${FSLDIR}/bin/calc_grad_perc_dev --fullwarp=${datadir}/warped/fullWarp -o ${datadir}/grad_dev
     ${FSLDIR}/bin/fslmerge -t ${datadir}/grad_dev ${datadir}/grad_dev_x ${datadir}/grad_dev_y ${datadir}/grad_dev_z
     ${FSLDIR}/bin/fslmaths ${datadir}/grad_dev -div 100 ${datadir}/grad_dev #Convert from % deviation to absolute
     ${FSLDIR}/bin/imrm ${datadir}/grad_dev_?
     ${FSLDIR}/bin/imrm ${datadir}/trilinear
-    ${FSLDIR}/bin/imrm ${datadir}/data_warped_vol1
-    ${FSLDIR}/bin/imrm ${datadir}/data_warped_dilated
-
-    #Keep the original warped data and warp fields
-    mkdir -p ${datadir}/warped
-    ${FSLDIR}/bin/immv ${datadir}/data_warped ${datadir}/warped
-    ${FSLDIR}/bin/immv ${datadir}/fov_mask_warped ${datadir}/warped
-    ${FSLDIR}/bin/immv ${datadir}/fullWarp ${datadir}/warped
-    ${FSLDIR}/bin/immv ${datadir}/fullWarp_abs ${datadir}/warped
+    ${FSLDIR}/bin/imrm ${datadir}/warped/data_dilated_vol1
 fi
 
 # mask out any data outside the field of view
 ${FSLDIR}/bin/fslmaths ${datadir}/data -mas ${datadir}/fov_mask ${datadir}/data
+${FSLDIR}/bin/fslmaths ${datadir}/cnr_maps -mas ${datadir}/fov_mask ${datadir}/cnr_maps
 
 # Remove negative intensity values (from eddy) from final data
 ${FSLDIR}/bin/fslmaths ${datadir}/data -thr 0 ${datadir}/data

DiffusionPreprocessing/scripts/run_eddy.sh

@@ -576,6 +576,7 @@ main()
 	eddy_command+="${sep_offs_moveOption} "
 	eddy_command+="${rmsOption} "
 	eddy_command+="${ff_valOption} "
+	eddy_command+="--cnr_maps "
 	eddy_command+="--imain=${workingdir}/Pos_Neg "
 	eddy_command+="--mask=${workingdir}/nodif_brain_mask "
 	eddy_command+="--index=${workingdir}/index.txt "

DiffusionPreprocessing/scripts/run_topup.sh

@@ -10,7 +10,7 @@ configdir=${HCPPIPEDIR_Config}
 #topup_config_file=${FSLDIR}/etc/flirtsch/b02b0.cnf
 topup_config_file=${configdir}/b02b0.cnf
 
-${FSLDIR}/bin/topup --imain=${workingdir}/Pos_Neg_b0 --datain=${workingdir}/acqparams.txt --config=${topup_config_file} --out=${workingdir}/topup_Pos_Neg_b0 -v
+${FSLDIR}/bin/topup --imain=${workingdir}/Pos_Neg_b0 --datain=${workingdir}/acqparams.txt --config=${topup_config_file} --out=${workingdir}/topup_Pos_Neg_b0 -v --fout=${workingdir}/topup_Pos_Neg_b0_field.nii.gz
 
 dimt=`${FSLDIR}/bin/fslval ${workingdir}/Pos_b0 dim4`
 dimt=$((${dimt} + 1))