diff --git a/000_view_preprocessing_results b/000_view_preprocessing_results new file mode 100644 index 0000000..03c79d5 --- /dev/null +++ b/000_view_preprocessing_results @@ -0,0 +1,67 @@ +Check_preprocessing () { + +subj=$1 + +low_res=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/$subj +high_res=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj + +# View denoising maps. The less anatomical structures are visible the better +mrview $low_res/Temporal/Distortions/$subj\_series* + +# Compare DWI before and after denoising for each encoding direction +mrview $low_res/Temporal/*AP.mif $low_res/Temporal/*AP_denoised.mif +mrview $low_res/Temporal/*PA.mif $low_res/Temporal/*PA_denoised.mif + +# Compare mrtrix gradient for corrected DWI with original gradient table +mrinfo $low_res/Temporal/$subj\_*AP_denoised.mif -export_grad_fsl $low_res/Temporal/bvecs $low_res/Temporal/bvals +gedit $low_res/Temporal/bvecs & +gedit $low_res/Temporal/Distortions/$subj\_gradients & +gedit /home/auguser2016/dMRI_DATA/RAW_DATA/$subj*/s*/*ep2d*/*.bvec +rm $low_res/Temporal/bvecs $low_res/Temporal/bvals + +# Compare DWI before and after dwipreproc +mrview $low_res/Temporal/*AP.mif $low_res/Temporal/*PA.mif $low_res/Temporal/$subj\_denoised_preproc.mif + +# View mask used for bias field correction +mrview $low_res/Temporal/$subj\_denoised_preproc.mif -overlay.load $low_res/$subj\_clean_150mm_mask.mif -overlay.opacity 0.2 -overlay.colourmap 3 + +# View bias field +mrview $low_res/Temporal/Distortions/$subj\_bias_field.mif + +# View DWI before and after bias correction +mrview $low_res/Temporal/$subj\_denoised_preproc.mif $low_res/$subj\_clean_150mm.mif + +# View normal and upsampled images +mrview $low_res/$subj\_clean_150mm.mif $high_res/$subj\_clean_075mm.mif + +# View upsampled mask +mrview $high_res/$subj\_clean_075mm.mif -overlay.load $high_res/$subj\_clean_075mm_mask.mif -overlay.opacity 0.2 -overlay.colourmap 3 + +# View SFR estimation voxels +mrview $high_res/$subj\_clean_075mm.mif -overlay.load $high_res/Temporal/$subj\_clean_075mm_SFR_voxels.mif + +# View SFR +shview $high_res/$subj\_clean_075mm_SFR.txt + +# View anatomical image overlaid on upsampled FOD +mrview $high_res/$subj\_clean_075mm.mif $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz + +# View results of bet2 +mrview $high_res/$subj\_anatomical.nii.gz -overlay.load $high_res/Temporal/$subj\_anatomical_brain.nii.gz -overlay.opacity 0.3 + +# View segmentation +mrview $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.load $high_res/$subj\_075mm_5tt.nii.gz -overlay.opacity 0.3 + +# Correct WM segmentation in chiasm - modify opened ROI and close mrview. Modified WM segmentation image will be shown - repeat previous steps until satisfactory result is obtained. Terminate console to finish process +mrview $high_res/$subj\_075mm_5tt.nii.gz -colourmap 3 -overlay.load $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $high_res/$subj\_5tt_patch.mif + +5ttedit $high_res/$subj\_075mm_5tt.nii.gz -wm $high_res/$subj\_5tt_patch.mif $high_res/$subj\_5tt_for_act_modified.nii.gz -force + +while true; do + mrview $high_res/$subj\_5tt_for_act_modified.nii.gz -colourmap 3 -overlay.load $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $high_res/$subj\_5tt_patch.mif + 5ttedit $high_res/$subj\_075mm_5tt.nii.gz -wm $high_res/$subj\_5tt_patch.mif $high_res/$subj\_5tt_for_act_modified.nii.gz -force +done + +} + +for i in /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/*; do Check_preprocessing ${i: -4}; done diff --git a/00_fsl_preprocess.sh b/00_fsl_preprocess.sh new file mode 100644 index 0000000..2f689b4 --- /dev/null +++ b/00_fsl_preprocess.sh @@ -0,0 +1,117 @@ +# Initializing + +# Folder with raw and preprocessed data +data=/home/auguser2016/dMRI_DATA + +Preparation () { + + # 2 resolutions (1.5 mm and 0.75mm), inside each folder all preprocessed data sets of corresponding resolution with masks, map of distortions etc. + a=${1::-6} + + low_res=PREPROCESSED_DATA/1.5mm_iso/$a + high_res=PREPROCESSED_DATA/0.75mm_iso/$a + + mkdir $data/$low_res + mkdir $data/$low_res/Temporal + mkdir $data/$low_res/Temporal/Distortions + + mkdir $data/$high_res + mkdir $data/$high_res/Temporal + + + cd $data/RAW_DATA/${1::-1}/s*/ + + # Raw data files - DWI and anatomy - unpacked with MRIcron (dcm2nii) as a 4D Nifti images. Raw data goes: RAW_DATA/data_set_name/study.../studies_X/ where important are those with ep2d + # Transforming Nifti into .mif files and denoising them. Denoising is the very first that should be done - it's results can be viewed. + + for b in *ep2d*; do + cd $b + mrconvert *a001.nii -stride -1,2,3,4 -fslgrad *.bvec *.bval $data/$low_res/Temporal/$a\_$b\.mif + dwidenoise $data/$low_res/Temporal/$a\_$b.mif $data/$low_res/Temporal/$a\_$b\_denoised.mif -noise $data/$low_res/Temporal/Distortions/$a\_$b\_noise.mif -debug + cd .. + done + + # Script from MRtrix interacting with fsl. Here version for correcting DWI with 2 phase encoding directions is used + dwipreproc AP $data/$low_res/Temporal/*AP_denoised.mif -rpe_all $data/$low_res/Temporal/*PA_denoised.mif $data/$low_res/Temporal/$a\_denoised_preproc.mif -export_grad_mrtrix $data/$low_res/Temporal/Distortions/$a\_gradients + + # Masking preprocessed data. Results should be checked + dwi2mask $data/$low_res/Temporal/$a\_denoised_preproc.mif $data/$low_res/$a\_clean_150mm_mask.mif + + # Masked is used for bias field correction using ANTS software + dwibiascorrect -ants $data/$low_res/Temporal/$a\_denoised_preproc.mif $data/$low_res/$a\_clean_150mm.mif -mask $data/$low_res/$a\_clean_150mm_mask.mif -bias $data/$low_res/Temporal/Distortions/$a\_bias_field.mif + + # Preprocessing of 1.5 mm iso data sets done. Moving to upsampled resolution + + # Upsampling by factor of 2 to 0.75mm. Justified step, in FBA even factor 3 can be recommended + mrresize $data/$low_res/$a\_clean_150mm.mif -scale 2.0 $data/$high_res/$a\_clean_075mm.mif + + # Calculating mask again + dwi2mask $data/$high_res/$a\_clean_075mm.mif $data/$high_res/$a\_clean_075mm_mask.mif + + # Estimating single fibre response (SFR) for default lmax + dwi2response tournier $data/$high_res/$a\_clean_075mm.mif $data/$high_res/$a\_clean_075mm_SFR.txt -shell 1600 -lmax 8 -mask $data/$high_res/$a\_clean_075mm_mask.mif -voxels $data/$high_res/Temporal/$a\_clean_075mm_SFR_voxels.mif + + # Estimating FOD for previously obtained SFR + dwiextract $data/$high_res/$a\_clean_075mm.mif - | dwi2fod msmt_csd - $data/$high_res/$a\_clean_075mm_SFR.txt $data/$high_res/$a\_clean_075mm_FOD.mif -mask $data/$high_res/$a\_clean_075mm_mask.mif + + # Performing coregistration of anatomical image to DWI images and segmenting it. It's easier approach as gradient table is not rotated this way. However for cross-modal studies it would be better to register DWI to anatomy + + cd *MPRAGE*iso + + # Converting DWI image from .mif to .nii.gz as FSL likes. This step should be checked as well. + mrconvert $data/$high_res/$a\_clean_075mm.mif $data/$high_res/Temporal/$a\_clean_075mm.nii.gz + + # Coping anatomical image to corresponding folder + cp o* $data/$high_res/$a\_anatomical.nii.gz + + # Extracting clean brain from anatomical image. Check this step + bet2 $data/$high_res/$a\_anatomical.nii.gz $data/$high_res/Temporal/$a\_anatomical_brain.nii.gz -f 0.6 + + # Corregister DWI to anatomy + epi_reg --epi=$data/$high_res/Temporal/$a\_clean_075mm.nii.gz --t1=$data/$high_res/$a\_anatomical.nii.gz --t1brain=$data/$high_res/Temporal/$a\_anatomical_brain.nii.gz --out=$data/$high_res/Temporal/$a\_nodif2mprage_2 --echospacing=0.00023 --pedir=-y + + # Reverse corregistration matrix + convert_xfm -inverse $data/$high_res/Temporal/$a\_nodif2mprage_2.mat -omat $data/$high_res/Temporal/$a\_mprage2nodif_2.mat + + # Apply reverted matrix to corregister anatomy to DWI + flirt -in $data/$high_res/$a\_anatomical.nii.gz -ref $data/$high_res/Temporal/$a\_clean_075mm.nii.gz -out $data/$high_res/$a\_075mm_anatomical_coreg2nodif_undist_2 -applyxfm -init $data/$high_res/Temporal/$a\_mprage2nodif_2.mat -interp sinc -sincwidth 7 -sincwindow hanning + + # Perform segmentation into WM, GM, subcortical GM, CSF and possible pathological tissue using FIRST from FSL + 5ttgen fsl $data/$high_res/$a\_075mm_anatomical_coreg2nodif_undist_2.nii.gz $data/$high_res/$a\_075mm_5tt.nii.gz + +} + +cd $data/RAW_DATA + +# List of all data_sets used in studies +list=(fe21_1325/ hw91_0844/ ib57_0731/ kw99_0633/ lw37_0977/ nb30_1185/ ow93_0974/ ps94_1516/ rx88_1234/ sj22_1218/ ta14_1065/ tq63_1214/ xs62_1217/ xn78_1085/ uh47_1309/ uf97_1072/ ow93_0974/ ) +#list=(fe21_1325/) + +# This way 4 sets are processed at once - we can't process all at once due to RAM usage, however preprocessing one at a time is ineffective because e.g. eddy from dwipreproc is using only one core. There is a version of eddy operating on multicore processors, but it's not running on this machine + +for i in ${list[@]:0:3}; do Preparation $i & done +wait +for i in ${list[@]:3:6}; do Preparation $i & done +wait +for i in ${list[@]:6:9}; do Preparation $i & done +wait +for i in ${list[@]:9:12}; do Preparation $i & done +wait +for i in ${list[@]:12:15}; do Preparation $i & done +wait +for i in ${list[@]:15:18}; do Preparation $i & done +wait +for i in ${list[@]:18:21}; do Preparation $i & done +wait +for i in ${list[@]:21:24}; do Preparation $i & done +wait +for i in ${list[@]:24:27}; do Preparation $i & done +wait +for i in ${list[@]:27:30}; do Preparation $i & done +wait +for i in ${list[@]:30:33}; do Preparation $i & done +wait + +# Returning to original folder +cd /home/auguser2016/Scripts/ + diff --git a/00_fsl_preprocess.sh~ b/00_fsl_preprocess.sh~ new file mode 100644 index 0000000..1196ef9 --- /dev/null +++ b/00_fsl_preprocess.sh~ @@ -0,0 +1,117 @@ +# Initializing + +# Folder with raw and preprocessed data +data=/home/auguser2016/dMRI_DATA + +Genral_preprocessing () { + + # 2 resolutions (1.5 mm and 0.75mm), inside each folder all preprocessed data sets of corresponding resolution with masks, map of distortions etc. + a=${1::-6} + + low_res=PREPROCESSED_DATA/1.5mm_iso/$a + high_res=PREPROCESSED_DATA/0.75mm_iso/$a + + mkdir $data/$low_res + mkdir $data/$low_res/Temporal + mkdir $data/$low_res/Temporal/Distortions + + mkdir $data/$high_res + mkdir $data/$high_res/Temporal + + + cd $data/RAW_DATA/${1::-1}/s*/ + + # Raw data files - DWI and anatomy - unpacked with MRIcron (dcm2nii) as a 4D Nifti images. Raw data goes: RAW_DATA/data_set_name/study.../studies_X/ where important are those with ep2d + # Transforming Nifti into .mif files and denoising them. Denoising is the very first that should be done - it's results can be viewed. + + for b in *ep2d*; do + cd $b + mrconvert *a001.nii -stride -1,2,3,4 -fslgrad *.bvec *.bval $data/$low_res/Temporal/$a\_$b\.mif + dwidenoise $data/$low_res/Temporal/$a\_$b.mif $data/$low_res/Temporal/$a\_$b\_denoised.mif -noise $data/$low_res/Temporal/Distortions/$a\_$b\_noise.mif -debug + cd .. + done + + # Script from MRtrix interacting with fsl. Here version for correcting DWI with 2 phase encoding directions is used + dwipreproc AP $data/$low_res/Temporal/*AP_denoised.mif -rpe_all $data/$low_res/Temporal/*PA_denoised.mif $data/$low_res/Temporal/$a\_denoised_preproc.mif -export_grad_mrtrix $data/$low_res/Temporal/Distortions/$a\_gradients + + # Masking preprocessed data. Results should be checked + dwi2mask $data/$low_res/Temporal/$a\_denoised_preproc.mif $data/$low_res/$a\_clean_150mm_mask.mif + + # Masked is used for bias field correction using ANTS software + dwibiascorrect -ants $data/$low_res/Temporal/$a\_denoised_preproc.mif $data/$low_res/$a\_clean_150mm.mif -mask $data/$low_res/$a\_clean_150mm_mask.mif -bias $data/$low_res/Temporal/Distortions/$a\_bias_field.mif + + # Preprocessing of 1.5 mm iso data sets done. Moving to upsampled resolution + + # Upsampling by factor of 2 to 0.75mm. Justified step, in FBA even factor 3 can be recommended + mrresize $data/$low_res/$a\_clean_150mm.mif -scale 2.0 $data/$high_res/$a\_clean_075mm.mif + + # Calculating mask again + dwi2mask $data/$high_res/$a\_clean_075mm.mif $data/$high_res/$a\_clean_075mm_mask.mif + + # Estimating single fibre response (SFR) for default lmax + dwi2response tournier $data/$high_res/$a\_clean_075mm.mif $data/$high_res/$a\_clean_075mm_SFR.txt -shell 1600 -lmax 8 -mask $data/$high_res/$a\_clean_075mm_mask.mif -voxels $data/$high_res/Temporal/$a\_clean_075mm_SFR_voxels.mif + + # Estimating FOD for previously obtained SFR + dwiextract $data/$high_res/$a\_clean_075mm.mif - | dwi2fod msmt_csd - $data/$high_res/$a\_clean_075mm_SFR.txt $data/$high_res/$a\_clean_075mm_FOD.mif -mask $data/$high_res/$a\_clean_075mm_mask.mif + + # Performing coregistration of anatomical image to DWI images and segmenting it. It's easier approach as gradient table is not rotated this way. However for cross-modal studies it would be better to register DWI to anatomy + + cd *MPRAGE*iso + + # Converting DWI image from .mif to .nii.gz as FSL likes. This step should be checked as well. + mrconvert $data/$high_res/$a\_clean_075mm.mif $data/$high_res/Temporal/$a\_clean_075mm.nii.gz + + # Coping anatomical image to corresponding folder + cp o* $data/$high_res/$a\_anatomical.nii.gz + + # Extracting clean brain from anatomical image. Check this step + bet2 $data/$high_res/$a\_anatomical.nii.gz $data/$high_res/Temporal/$a\_anatomical_brain.nii.gz -f 0.6 + + # Corregister DWI to anatomy + epi_reg --epi=$data/$high_res/Temporal/$a\_clean_075mm.nii.gz --t1=$data/$high_res/$a\_anatomical.nii.gz --t1brain=$data/$high_res/Temporal/$a\_anatomical_brain.nii.gz --out=$data/$high_res/Temporal/$a\_nodif2mprage_2 --echospacing=0.00023 --pedir=-y + + # Reverse corregistration matrix + convert_xfm -inverse $data/$high_res/Temporal/$a\_nodif2mprage_2.mat -omat $data/$high_res/Temporal/$a\_mprage2nodif_2.mat + + # Apply reverted matrix to corregister anatomy to DWI + flirt -in $data/$high_res/$a\_anatomical.nii.gz -ref $data/$high_res/Temporal/$a\_clean_075mm.nii.gz -out $data/$high_res/$a\_075mm_anatomical_coreg2nodif_undist_2 -applyxfm -init $data/$high_res/Temporal/$a\_mprage2nodif_2.mat -interp sinc -sincwidth 7 -sincwindow hanning + + # Perform segmentation into WM, GM, subcortical GM, CSF and possible pathological tissue using FIRST from FSL + 5ttgen fsl $data/$high_res/$a\_075mm_anatomical_coreg2nodif_undist_2.nii.gz $data/$high_res/$a\_075mm_5tt.nii.gz + +} + +cd $data/RAW_DATA + +# List of all data_sets used in studies +list=(fe21_1325/ hw91_0844/ ib57_0731/ kw99_0633/ lw37_0977/ nb30_1185/ ow93_0974/ ps94_1516/ rx88_1234/ sj22_1218/ ta14_1065/ tq63_1214/ xs62_1217/ xn78_1085/ uh47_1309/ uf97_1072/ ow93_0974/ ) +#list=(fe21_1325/) + +# This way 4 sets are processed at once - we can't process all at once due to RAM usage, however preprocessing one at a time is ineffective because e.g. eddy from dwipreproc is using only one core. There is a version of eddy operating on multicore processors, but it's not running on this machine + +for i in ${list[@]:0:3}; do Preparation $i & done +wait +for i in ${list[@]:3:6}; do Preparation $i & done +wait +for i in ${list[@]:6:9}; do Preparation $i & done +wait +for i in ${list[@]:9:12}; do Preparation $i & done +wait +for i in ${list[@]:12:15}; do Preparation $i & done +wait +for i in ${list[@]:15:18}; do Preparation $i & done +wait +for i in ${list[@]:18:21}; do Preparation $i & done +wait +for i in ${list[@]:21:24}; do Preparation $i & done +wait +for i in ${list[@]:24:27}; do Preparation $i & done +wait +for i in ${list[@]:27:30}; do Preparation $i & done +wait +for i in ${list[@]:30:33}; do Preparation $i & done +wait + +# Returning to original folder +cd /home/auguser2016/Scripts/ + diff --git a/00_preprocess.sh b/00_preprocess.sh new file mode 100644 index 0000000..67b2603 --- /dev/null +++ b/00_preprocess.sh @@ -0,0 +1,118 @@ +# Initializing + +# Folder with raw and preprocessed data +data=/home/auguser2016/dMRI_DATA + +General_preprocessing () { + + # 2 resolutions (1.5 mm and 0.75mm), inside each folder all preprocessed data sets of corresponding resolution with masks, map of distortions etc. + a=${1::-6} + + low_res=PREPROCESSED_DATA/1.5mm_iso/$a + high_res=PREPROCESSED_DATA/0.75mm_iso/$a + + mkdir $data/$low_res + mkdir $data/$low_res/Temporal + mkdir $data/$low_res/Temporal/Distortions + + mkdir $data/$high_res + mkdir $data/$high_res/Temporal + + + cd $data/RAW_DATA/${1::-1}/s*/ + + # Raw data files - DWI and anatomy - unpacked with MRIcron (dcm2nii) as a 4D Nifti images. Raw data goes: RAW_DATA/data_set_name/study.../studies_X/ where important are those with ep2d + # Transforming Nifti into .mif files and denoising them. Denoising is the very first that should be done - it's results can be viewed. + + for b in *ep2d*; do + cd $b + mrconvert *a001.nii -stride -1,2,3,4 -fslgrad *.bvec *.bval $data/$low_res/Temporal/$a\_$b\.mif + dwidenoise $data/$low_res/Temporal/$a\_$b.mif $data/$low_res/Temporal/$a\_$b\_denoised.mif -noise $data/$low_res/Temporal/Distortions/$a\_$b\_noise.mif -debug + cd .. + done + + # Script from MRtrix interacting with fsl. Here version for correcting DWI with 2 phase encoding directions is used + dwipreproc AP $data/$low_res/Temporal/*AP_denoised.mif -rpe_all $data/$low_res/Temporal/*PA_denoised.mif $data/$low_res/Temporal/$a\_denoised_preproc.mif -export_grad_mrtrix $data/$low_res/Temporal/Distortions/$a\_gradients + + # Masking preprocessed data. Results should be checked + dwi2mask $data/$low_res/Temporal/$a\_denoised_preproc.mif $data/$low_res/$a\_clean_150mm_mask.mif + + # Masked is used for bias field correction using ANTS software + dwibiascorrect -ants $data/$low_res/Temporal/$a\_denoised_preproc.mif $data/$low_res/$a\_clean_150mm.mif -mask $data/$low_res/$a\_clean_150mm_mask.mif -bias $data/$low_res/Temporal/Distortions/$a\_bias_field.mif + + # Preprocessing of 1.5 mm iso data sets done. Moving to upsampled resolution + + # Upsampling by factor of 2 to 0.75mm. Justified step, in FBA even factor 3 can be recommended + mrresize $data/$low_res/$a\_clean_150mm.mif -scale 2.0 $data/$high_res/$a\_clean_075mm.mif + + # Calculating mask again + dwi2mask $data/$high_res/$a\_clean_075mm.mif $data/$high_res/$a\_clean_075mm_mask.mif + + # Estimating single fibre response (SFR) for default lmax + dwi2response tournier $data/$high_res/$a\_clean_075mm.mif $data/$high_res/$a\_clean_075mm_SFR.txt -shell 1600 -lmax 8 -mask $data/$high_res/$a\_clean_075mm_mask.mif -voxels $data/$high_res/Temporal/$a\_clean_075mm_SFR_voxels.mif + + # Estimating FOD for previously obtained SFR + dwiextract $data/$high_res/$a\_clean_075mm.mif - | dwi2fod msmt_csd - $data/$high_res/$a\_clean_075mm_SFR.txt $data/$high_res/$a\_clean_075mm_FOD.mif -mask $data/$high_res/$a\_clean_075mm_mask.mif + + # Performing coregistration of anatomical image to DWI images and segmenting it. It's easier approach as gradient table is not rotated this way. However for cross-modal studies it would be better to register DWI to anatomy + + cd *MPRAGE*iso + + # Converting DWI image from .mif to .nii.gz as FSL likes. This step should be checked as well. + mrconvert $data/$high_res/$a\_clean_075mm.mif $data/$high_res/Temporal/$a\_clean_075mm.nii.gz + + # Coping anatomical image to corresponding folder + cp o* $data/$high_res/$a\_anatomical.nii.gz + + # Extracting clean brain from anatomical image. Check this step + bet2 $data/$high_res/$a\_anatomical.nii.gz $data/$high_res/Temporal/$a\_anatomical_brain.nii.gz -f 0.6 + + # Corregister DWI to anatomy + epi_reg --epi=$data/$high_res/Temporal/$a\_clean_075mm.nii.gz --t1=$data/$high_res/$a\_anatomical.nii.gz --t1brain=$data/$high_res/Temporal/$a\_anatomical_brain.nii.gz --out=$data/$high_res/Temporal/$a\_nodif2mprage_2 --echospacing=0.00023 --pedir=-y + + # Reverse corregistration matrix + convert_xfm -inverse $data/$high_res/Temporal/$a\_nodif2mprage_2.mat -omat $data/$high_res/Temporal/$a\_mprage2nodif_2.mat + + # Apply reverted matrix to corregister anatomy to DWI + flirt -in $data/$high_res/$a\_anatomical.nii.gz -ref $data/$high_res/Temporal/$a\_clean_075mm.nii.gz -out $data/$high_res/$a\_075mm_anatomical_coreg2nodif_undist_2 -applyxfm -init $data/$high_res/Temporal/$a\_mprage2nodif_2.mat -interp sinc -sincwidth 7 -sincwindow hanning + + # Perform segmentation into WM, GM, subcortical GM, CSF and possible pathological tissue using FIRST from FSL + 5ttgen fsl $data/$high_res/$a\_075mm_anatomical_coreg2nodif_undist_2.nii.gz $data/$high_res/$a\_075mm_5tt.nii.gz + +} + +cd $data/RAW_DATA + +# List of all data_sets used in studies +list=(fe21_1325/ hw91_0844/ ib57_0731/ kw99_0633/ lw37_0977/ nb30_1185/ ow93_0974/ ps94_1516/ rx88_1234/ sj22_1218/ ta14_1065/ tq63_1214/ xs62_1217/ xn78_1085/ uh47_1309/ uf97_1072/ ow93_0974/ ) + +list=(uh47_1309/) + +# This way 4 sets are processed at once - we can't process all at once due to RAM usage, however preprocessing one at a time is ineffective because e.g. eddy from dwipreproc is using only one core. There is a version of eddy operating on multicore processors, but it's not running on this machine + +for i in ${list[@]:0:3}; do General_preprocessing $i & done +wait +for i in ${list[@]:3:6}; do General_preprocessing $i & done +wait +for i in ${list[@]:6:9}; do General_preprocessing $i & done +wait +for i in ${list[@]:9:12}; do General_preprocessing $i & done +wait +for i in ${list[@]:12:15}; do General_preprocessing $i & done +wait +for i in ${list[@]:15:18}; do General_preprocessing $i & done +wait +for i in ${list[@]:18:21}; do General_preprocessing $i & done +wait +for i in ${list[@]:21:24}; do General_preprocessing $i & done +wait +for i in ${list[@]:24:27}; do General_preprocessing $i & done +wait +for i in ${list[@]:27:30}; do General_preprocessing $i & done +wait +for i in ${list[@]:30:33}; do General_preprocessing $i & done +wait + +# Returning to original folder +cd /home/auguser2016/Scripts/ + diff --git a/00_preprocess.sh~ b/00_preprocess.sh~ new file mode 100644 index 0000000..3985a09 --- /dev/null +++ b/00_preprocess.sh~ @@ -0,0 +1,117 @@ +# Initializing + +# Folder with raw and preprocessed data +data=/home/auguser2016/dMRI_DATA + +General_preprocessing () { + + # 2 resolutions (1.5 mm and 0.75mm), inside each folder all preprocessed data sets of corresponding resolution with masks, map of distortions etc. + a=${1::-6} + + low_res=PREPROCESSED_DATA/1.5mm_iso/$a + high_res=PREPROCESSED_DATA/0.75mm_iso/$a + + mkdir $data/$low_res + mkdir $data/$low_res/Temporal + mkdir $data/$low_res/Temporal/Distortions + + mkdir $data/$high_res + mkdir $data/$high_res/Temporal + + + cd $data/RAW_DATA/${1::-1}/s*/ + + # Raw data files - DWI and anatomy - unpacked with MRIcron (dcm2nii) as a 4D Nifti images. Raw data goes: RAW_DATA/data_set_name/study.../studies_X/ where important are those with ep2d + # Transforming Nifti into .mif files and denoising them. Denoising is the very first that should be done - it's results can be viewed. + + for b in *ep2d*; do + cd $b + mrconvert *a001.nii -stride -1,2,3,4 -fslgrad *.bvec *.bval $data/$low_res/Temporal/$a\_$b\.mif + dwidenoise $data/$low_res/Temporal/$a\_$b.mif $data/$low_res/Temporal/$a\_$b\_denoised.mif -noise $data/$low_res/Temporal/Distortions/$a\_$b\_noise.mif -debug + cd .. + done + + # Script from MRtrix interacting with fsl. Here version for correcting DWI with 2 phase encoding directions is used + dwipreproc AP $data/$low_res/Temporal/*AP_denoised.mif -rpe_all $data/$low_res/Temporal/*PA_denoised.mif $data/$low_res/Temporal/$a\_denoised_preproc.mif -export_grad_mrtrix $data/$low_res/Temporal/Distortions/$a\_gradients + + # Masking preprocessed data. Results should be checked + dwi2mask $data/$low_res/Temporal/$a\_denoised_preproc.mif $data/$low_res/$a\_clean_150mm_mask.mif + + # Masked is used for bias field correction using ANTS software + dwibiascorrect -ants $data/$low_res/Temporal/$a\_denoised_preproc.mif $data/$low_res/$a\_clean_150mm.mif -mask $data/$low_res/$a\_clean_150mm_mask.mif -bias $data/$low_res/Temporal/Distortions/$a\_bias_field.mif + + # Preprocessing of 1.5 mm iso data sets done. Moving to upsampled resolution + + # Upsampling by factor of 2 to 0.75mm. Justified step, in FBA even factor 3 can be recommended + mrresize $data/$low_res/$a\_clean_150mm.mif -scale 2.0 $data/$high_res/$a\_clean_075mm.mif + + # Calculating mask again + dwi2mask $data/$high_res/$a\_clean_075mm.mif $data/$high_res/$a\_clean_075mm_mask.mif + + # Estimating single fibre response (SFR) for default lmax + dwi2response tournier $data/$high_res/$a\_clean_075mm.mif $data/$high_res/$a\_clean_075mm_SFR.txt -shell 1600 -lmax 8 -mask $data/$high_res/$a\_clean_075mm_mask.mif -voxels $data/$high_res/Temporal/$a\_clean_075mm_SFR_voxels.mif + + # Estimating FOD for previously obtained SFR + dwiextract $data/$high_res/$a\_clean_075mm.mif - | dwi2fod msmt_csd - $data/$high_res/$a\_clean_075mm_SFR.txt $data/$high_res/$a\_clean_075mm_FOD.mif -mask $data/$high_res/$a\_clean_075mm_mask.mif + + # Performing coregistration of anatomical image to DWI images and segmenting it. It's easier approach as gradient table is not rotated this way. However for cross-modal studies it would be better to register DWI to anatomy + + cd *MPRAGE*iso + + # Converting DWI image from .mif to .nii.gz as FSL likes. This step should be checked as well. + mrconvert $data/$high_res/$a\_clean_075mm.mif $data/$high_res/Temporal/$a\_clean_075mm.nii.gz + + # Coping anatomical image to corresponding folder + cp o* $data/$high_res/$a\_anatomical.nii.gz + + # Extracting clean brain from anatomical image. Check this step + bet2 $data/$high_res/$a\_anatomical.nii.gz $data/$high_res/Temporal/$a\_anatomical_brain.nii.gz -f 0.6 + + # Corregister DWI to anatomy + epi_reg --epi=$data/$high_res/Temporal/$a\_clean_075mm.nii.gz --t1=$data/$high_res/$a\_anatomical.nii.gz --t1brain=$data/$high_res/Temporal/$a\_anatomical_brain.nii.gz --out=$data/$high_res/Temporal/$a\_nodif2mprage_2 --echospacing=0.00023 --pedir=-y + + # Reverse corregistration matrix + convert_xfm -inverse $data/$high_res/Temporal/$a\_nodif2mprage_2.mat -omat $data/$high_res/Temporal/$a\_mprage2nodif_2.mat + + # Apply reverted matrix to corregister anatomy to DWI + flirt -in $data/$high_res/$a\_anatomical.nii.gz -ref $data/$high_res/Temporal/$a\_clean_075mm.nii.gz -out $data/$high_res/$a\_075mm_anatomical_coreg2nodif_undist_2 -applyxfm -init $data/$high_res/Temporal/$a\_mprage2nodif_2.mat -interp sinc -sincwidth 7 -sincwindow hanning + + # Perform segmentation into WM, GM, subcortical GM, CSF and possible pathological tissue using FIRST from FSL + 5ttgen fsl $data/$high_res/$a\_075mm_anatomical_coreg2nodif_undist_2.nii.gz $data/$high_res/$a\_075mm_5tt.nii.gz + +} + +cd $data/RAW_DATA + +# List of all data_sets used in studies +list=(fe21_1325/ hw91_0844/ ib57_0731/ kw99_0633/ lw37_0977/ nb30_1185/ ow93_0974/ ps94_1516/ rx88_1234/ sj22_1218/ ta14_1065/ tq63_1214/ xs62_1217/ xn78_1085/ uh47_1309/ uf97_1072/ ow93_0974/ ) + + +# This way 4 sets are processed at once - we can't process all at once due to RAM usage, however preprocessing one at a time is ineffective because e.g. eddy from dwipreproc is using only one core. There is a version of eddy operating on multicore processors, but it's not running on this machine + +for i in ${list[@]:0:3}; do General_preprocessing $i & done +wait +for i in ${list[@]:3:6}; do General_preprocessing $i & done +wait +for i in ${list[@]:6:9}; do General_preprocessing $i & done +wait +for i in ${list[@]:9:12}; do General_preprocessing $i & done +wait +for i in ${list[@]:12:15}; do General_preprocessing $i & done +wait +for i in ${list[@]:15:18}; do General_preprocessing $i & done +wait +for i in ${list[@]:18:21}; do General_preprocessing $i & done +wait +for i in ${list[@]:21:24}; do General_preprocessing $i & done +wait +for i in ${list[@]:24:27}; do General_preprocessing $i & done +wait +for i in ${list[@]:27:30}; do General_preprocessing $i & done +wait +for i in ${list[@]:30:33}; do General_preprocessing $i & done +wait + +# Returning to original folder +cd /home/auguser2016/Scripts/ + diff --git a/01_check_preprocessing.sh b/01_check_preprocessing.sh new file mode 100644 index 0000000..8b717c0 --- /dev/null +++ b/01_check_preprocessing.sh @@ -0,0 +1,81 @@ +Check_preprocessing () { + +subj=$1 + +low_res=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/$subj +high_res=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj + +# View denoising maps. The less anatomical structures are visible the better +mrview $low_res/Temporal/Distortions/$subj\_series* + +# Compare DWI before and after denoising for each encoding direction +mrview $low_res/Temporal/*AP.mif $low_res/Temporal/*AP_denoised.mif +mrview $low_res/Temporal/*PA.mif $low_res/Temporal/*PA_denoised.mif + +# Compare mrtrix gradient for corrected DWI with original gradient table +mrinfo $low_res/Temporal/$subj\_*AP_denoised.mif -export_grad_fsl $low_res/Temporal/bvecs $low_res/Temporal/bvals +gedit $low_res/Temporal/bvecs & +gedit $low_res/Temporal/Distortions/$subj\_gradients & +gedit /home/auguser2016/dMRI_DATA/RAW_DATA/$subj*/s*/*ep2d*/*.bvec +rm $low_res/Temporal/bvecs $low_res/Temporal/bvals + +# Compare DWI before and after dwipreproc +mrview $low_res/Temporal/*AP.mif $low_res/Temporal/*PA.mif $low_res/Temporal/$subj\_denoised_preproc.mif + +# View mask used for bias field correction +mrview $low_res/Temporal/$subj\_denoised_preproc.mif -overlay.load $low_res/$subj\_clean_150mm_mask.mif -overlay.opacity 0.2 -overlay.colourmap 3 + +# View bias field +mrview $low_res/Temporal/Distortions/$subj\_bias_field.mif + +# View DWI before and after bias correction +mrview $low_res/Temporal/$subj\_denoised_preproc.mif $low_res/$subj\_clean_150mm.mif + +# View normal and upsampled images +mrview $low_res/$subj\_clean_150mm.mif $high_res/$subj\_clean_075mm.mif + +# View upsampled mask +mrview $high_res/$subj\_clean_075mm.mif -overlay.load $high_res/$subj\_clean_075mm_mask.mif -overlay.opacity 0.2 -overlay.colourmap 3 + +# View SFR estimation voxels +mrview $high_res/$subj\_clean_075mm.mif -overlay.load $high_res/Temporal/$subj\_clean_075mm_SFR_voxels.mif + +# View SFR +shview $high_res/$subj\_clean_075mm_SFR.txt + +# View anatomical image overlaid on DWI +mrview $high_res/$subj\_clean_075mm_FOD.mif $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz + +# View results of bet2 +mrview $high_res/$subj\_anatomical.nii.gz -overlay.load $high_res/Temporal/$subj\_anatomical_brain.nii.gz -overlay.opacity 0.3 + +# View segmentation +mrview $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.load $high_res/$subj\_075mm_5tt.nii.gz -overlay.opacity 0.3 + +# Correct WM segmentation in chiasm - modify opened ROI and close mrview. Modified WM segmentation image will be shown - repeat previous steps until satisfactory result is obtained. Terminate console to finish process +mrview $high_res/$subj\_075mm_5tt.nii.gz -colourmap 3 -overlay.load $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $high_res/$subj\_5tt_patch.mif + +5ttedit $high_res/$subj\_075mm_5tt.nii.gz -wm $high_res/$subj\_5tt_patch.mif $high_res/$subj\_5tt_for_act_modified.nii.gz -force + +while true; do + mrview $high_res/$subj\_5tt_for_act_modified.nii.gz -colourmap 3 -overlay.load $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $high_res/$subj\_5tt_patch.mif + 5ttedit $high_res/$subj\_075mm_5tt.nii.gz -wm $high_res/$subj\_5tt_patch.mif $high_res/$subj\_5tt_for_act_modified.nii.gz -force +done + +#data=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj + +### Correct WM segmentation in chiasm - modify opened ROI and close mrview. Modified WM segmentation image will be shown - repeat previous steps until satisfactory result is obtained. Terminate console to finish process +#mrview $data/$subj\_075mm_5tt.nii.gz -colourmap 3 -overlay.load $data/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $data/$subj\_5tt_patch.mif + +#5ttedit $data/$subj\_075mm_5tt.nii.gz -wm $data/$subj\_5tt_patch.mif $data/$subj\_075mm_5tt_modify.nii.gz -force + +#while true; do +# mrview $data/$subj\_075mm_5tt_modify.nii.gz -colourmap 3 -overlay.load $data/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $data/$subj\_5tt_patch.mif +# 5ttedit $data/$subj\_075mm_5tt.nii.gz -wm $data/$subj\_5tt_patch.mif $data/$subj\_075mm_5tt_modify.nii.gz -force +#done + +} + +#for i in /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/*; do Check_preprocessing ${i: -4}; done + +Check_preprocessing $1 diff --git a/01_check_preprocessing.sh~ b/01_check_preprocessing.sh~ new file mode 100644 index 0000000..8b717c0 --- /dev/null +++ b/01_check_preprocessing.sh~ @@ -0,0 +1,81 @@ +Check_preprocessing () { + +subj=$1 + +low_res=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/$subj +high_res=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj + +# View denoising maps. The less anatomical structures are visible the better +mrview $low_res/Temporal/Distortions/$subj\_series* + +# Compare DWI before and after denoising for each encoding direction +mrview $low_res/Temporal/*AP.mif $low_res/Temporal/*AP_denoised.mif +mrview $low_res/Temporal/*PA.mif $low_res/Temporal/*PA_denoised.mif + +# Compare mrtrix gradient for corrected DWI with original gradient table +mrinfo $low_res/Temporal/$subj\_*AP_denoised.mif -export_grad_fsl $low_res/Temporal/bvecs $low_res/Temporal/bvals +gedit $low_res/Temporal/bvecs & +gedit $low_res/Temporal/Distortions/$subj\_gradients & +gedit /home/auguser2016/dMRI_DATA/RAW_DATA/$subj*/s*/*ep2d*/*.bvec +rm $low_res/Temporal/bvecs $low_res/Temporal/bvals + +# Compare DWI before and after dwipreproc +mrview $low_res/Temporal/*AP.mif $low_res/Temporal/*PA.mif $low_res/Temporal/$subj\_denoised_preproc.mif + +# View mask used for bias field correction +mrview $low_res/Temporal/$subj\_denoised_preproc.mif -overlay.load $low_res/$subj\_clean_150mm_mask.mif -overlay.opacity 0.2 -overlay.colourmap 3 + +# View bias field +mrview $low_res/Temporal/Distortions/$subj\_bias_field.mif + +# View DWI before and after bias correction +mrview $low_res/Temporal/$subj\_denoised_preproc.mif $low_res/$subj\_clean_150mm.mif + +# View normal and upsampled images +mrview $low_res/$subj\_clean_150mm.mif $high_res/$subj\_clean_075mm.mif + +# View upsampled mask +mrview $high_res/$subj\_clean_075mm.mif -overlay.load $high_res/$subj\_clean_075mm_mask.mif -overlay.opacity 0.2 -overlay.colourmap 3 + +# View SFR estimation voxels +mrview $high_res/$subj\_clean_075mm.mif -overlay.load $high_res/Temporal/$subj\_clean_075mm_SFR_voxels.mif + +# View SFR +shview $high_res/$subj\_clean_075mm_SFR.txt + +# View anatomical image overlaid on DWI +mrview $high_res/$subj\_clean_075mm_FOD.mif $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz + +# View results of bet2 +mrview $high_res/$subj\_anatomical.nii.gz -overlay.load $high_res/Temporal/$subj\_anatomical_brain.nii.gz -overlay.opacity 0.3 + +# View segmentation +mrview $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.load $high_res/$subj\_075mm_5tt.nii.gz -overlay.opacity 0.3 + +# Correct WM segmentation in chiasm - modify opened ROI and close mrview. Modified WM segmentation image will be shown - repeat previous steps until satisfactory result is obtained. Terminate console to finish process +mrview $high_res/$subj\_075mm_5tt.nii.gz -colourmap 3 -overlay.load $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $high_res/$subj\_5tt_patch.mif + +5ttedit $high_res/$subj\_075mm_5tt.nii.gz -wm $high_res/$subj\_5tt_patch.mif $high_res/$subj\_5tt_for_act_modified.nii.gz -force + +while true; do + mrview $high_res/$subj\_5tt_for_act_modified.nii.gz -colourmap 3 -overlay.load $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $high_res/$subj\_5tt_patch.mif + 5ttedit $high_res/$subj\_075mm_5tt.nii.gz -wm $high_res/$subj\_5tt_patch.mif $high_res/$subj\_5tt_for_act_modified.nii.gz -force +done + +#data=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj + +### Correct WM segmentation in chiasm - modify opened ROI and close mrview. Modified WM segmentation image will be shown - repeat previous steps until satisfactory result is obtained. Terminate console to finish process +#mrview $data/$subj\_075mm_5tt.nii.gz -colourmap 3 -overlay.load $data/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $data/$subj\_5tt_patch.mif + +#5ttedit $data/$subj\_075mm_5tt.nii.gz -wm $data/$subj\_5tt_patch.mif $data/$subj\_075mm_5tt_modify.nii.gz -force + +#while true; do +# mrview $data/$subj\_075mm_5tt_modify.nii.gz -colourmap 3 -overlay.load $data/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $data/$subj\_5tt_patch.mif +# 5ttedit $data/$subj\_075mm_5tt.nii.gz -wm $data/$subj\_5tt_patch.mif $data/$subj\_075mm_5tt_modify.nii.gz -force +#done + +} + +#for i in /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/*; do Check_preprocessing ${i: -4}; done + +Check_preprocessing $1 diff --git a/01_preparation_for_acpc.sh~ b/01_preparation_for_acpc.sh~ new file mode 100644 index 0000000..dba2704 --- /dev/null +++ b/01_preparation_for_acpc.sh~ @@ -0,0 +1,22 @@ +ACPC_preparation () { + + # Registration to ACPC space requires preprocessed DWI image in .nii (or .nii.gz) format with corresponding bvecs and bvals. T1 image in ACPC space will be prepared separately + + # Defining data set and its' path + a=${1::-6} + data=/home/auguser2016/dMRI_DATA + high_res=PREPROCESSED_DATA/0.75mm_iso/$a + + # Obtain bvecs and bvals + mrinfo $data/$high_res/$a\_clean_075mm.mif -export_grad_fsl $data/$high_res/$a\_bvecs.bvecs $data/$high_res/$a\_bvals.bvals + + # Transform DW in .mif into .nii.gz + mrconvert $data/$high_res/$a\_clean_075mm.mif $data/$high_res/$a\_clean_075mm.nii.gz +} + +# List of all data_sets used in studies +list=(fe21_1325/ hw91_0844/ ib57_0731/ kw99_0633/ lw37_0977/ nb30_1185/ ow93_0974/ ps94_1516/ rx88_1234/ sj22_1218/ ta14_1065/ tq63_1214/ xs62_1217/ xn78_1085/ uh47_1309/ uf97_1072/ ow93_0974/ ) +list=(fe21_1325/) + +# As during processing of each data set full resources are used, lining them up is an effective approach +for i in ${list[@]}; do ACPC_preparation $i; done diff --git a/01_prepare_for_acpc.sh b/01_prepare_for_acpc.sh new file mode 100644 index 0000000..ca1d232 --- /dev/null +++ b/01_prepare_for_acpc.sh @@ -0,0 +1,21 @@ +ACPC_preparation () { + + # Registration to ACPC space requires preprocessed DWI image in .nii (or .nii.gz) format with corresponding bvecs and bvals. T1 image in ACPC space will be prepared separately + + # Defining data set and its' path + a=${1::-6} + data=/home/auguser2016/dMRI_DATA + high_res=PREPROCESSED_DATA/0.75mm_iso/$a + + # Obtain bvecs and bvals + mrinfo $data/$high_res/$a\_clean_075mm.mif -export_grad_fsl $data/$high_res/$a\_bvecs.bvecs $data/$high_res/$a\_bvals.bvals + + # Transform DW in .mif into .nii.gz + mrconvert $data/$high_res/$a\_clean_075mm.mif $data/$high_res/$a\_clean_075mm.nii.gz +} + +# List of all data_sets used in studies +list=(fe21_1325/ hw91_0844/ ib57_0731/ kw99_0633/ lw37_0977/ nb30_1185/ ow93_0974/ ps94_1516/ rx88_1234/ sj22_1218/ ta14_1065/ tq63_1214/ xs62_1217/ xn78_1085/ uh47_1309/ uf97_1072/ ow93_0974/ ) + +# As during processing of each data set full resources are used, lining them up is an effective approach +for i in ${list[@]}; do ACPC_preparation $i; done diff --git a/testing_acpc.m b/01_testing_acpc.m similarity index 100% rename from testing_acpc.m rename to 01_testing_acpc.m diff --git a/02_ensemble_tracking.txt b/02_ensemble_tracking.txt new file mode 100644 index 0000000..c5223eb --- /dev/null +++ b/02_ensemble_tracking.txt @@ -0,0 +1,63 @@ +# 1 Argument: ID of participant + +# Requirements: +# Data preprocessed +# DWI corregisterd to anatomy +# Masks estimated +# ROIs drawn + +subj=$1 +analysis_folder=/media/auguser2016/Volume/Quantifying_Chiasm + +mkdir $analysis_folder/$subj + +output=$analysis_folder/$subj + +# SFR Estimation +lmax=6 'This should be as low as possible, optimal for most participants value #SFR for lower lmax, keep it and use higher 8,10,12 for fod estimation OR Get SFR with optimal lmax for most participants and use given lmax for all participants. Lower lmax preffered' +dwi2response tournier $subj\_complete.mif -shell 1600 -lmax $lmax -mask $subj\_mask_upsampled.nii.gz $output/$subj\_sfr.txt + +# Estimation of Fiber Orientation Distribution function +for lmax in "8 10 12"; do + dwiextract $subj\_complete.mif - | dwi2fod msmt_csd - $output/$subj\_sfr.txt $output/$subj\_fod_lmax_$lmax\.nii.gz -mask $subj\_mask_upsampled.nii.gz -lmax $lmax +done + +# Generation of 5TT image based on acpc registerred image +5ttgen fsl anatomy/$subj\_t1_acpc.nii.gz $subj\_5tt.nii.gz + +# Generate tractogram + +lmax="8 10 12" +FAs="0.05" +Curv="30 45 60" + +for i in $lmax; do + for j in $FAs; do + for k in $Curv; do + tckgen $subj\_fod_lmax_$lmax\.nii.gz lw_results_r_r_$i\_$j\_$k.tck -number 300 -cutoff $j -angle $k -act lw_5tt.nii.gz -unidirectional -seed_image $subj\_start_right.mif -include $subj\_end_right.mif -stop -step 0.2 -rk4; + tckgen $subj\_fod_lmax_$lmax\.nii.gz lw_results_r_l_$i\_$j\_$k.tck -number 300 -cutoff $j -angle $k -act lw_5tt.nii.gz -unidirectional -seed_image $subj\_start_right.mif -include $subj\_end_left.mif -stop -step 0.2 -rk4; + tckgen $subj\_fod_lmax_$lmax\.nii.gz lw_results_l_r_$i\_$j\_$k.tck -number 300 -cutoff $j -angle $k -act lw_5tt.nii.gz -unidirectional -seed_image $subj\_start_left.mif -include $subj\_end_right.mif -stop -step 0.2 -rk4; + tckgen $subj\_fod_lmax_$lmax\.nii.gz lw_results_l_l_$i\_$j\_$k.tck -number 300 -cutoff $j -angle $k -act lw_5tt.nii.gz -unidirectional -seed_image $subj\_start_left.mif -include $subj\_end_left.mif -stop -step 0.2 -rk4; + echo " " + done + done +done + +altogether="" + +for z in *.tck; do + altogether+="${z} " +done + +tckedit $altogether $subj\_et.tck -force + + +source 03_life.m $subj + +#for r in $subjects; do +# analizuj "$r" +# done + +matlab -nodisplay -r 03_life.m + + diff --git a/02_ensemble_tracking.txt~ b/02_ensemble_tracking.txt~ new file mode 100644 index 0000000..c5223eb --- /dev/null +++ b/02_ensemble_tracking.txt~ @@ -0,0 +1,63 @@ +# 1 Argument: ID of participant + +# Requirements: +# Data preprocessed +# DWI corregisterd to anatomy +# Masks estimated +# ROIs drawn + +subj=$1 +analysis_folder=/media/auguser2016/Volume/Quantifying_Chiasm + +mkdir $analysis_folder/$subj + +output=$analysis_folder/$subj + +# SFR Estimation +lmax=6 'This should be as low as possible, optimal for most participants value #SFR for lower lmax, keep it and use higher 8,10,12 for fod estimation OR Get SFR with optimal lmax for most participants and use given lmax for all participants. Lower lmax preffered' +dwi2response tournier $subj\_complete.mif -shell 1600 -lmax $lmax -mask $subj\_mask_upsampled.nii.gz $output/$subj\_sfr.txt + +# Estimation of Fiber Orientation Distribution function +for lmax in "8 10 12"; do + dwiextract $subj\_complete.mif - | dwi2fod msmt_csd - $output/$subj\_sfr.txt $output/$subj\_fod_lmax_$lmax\.nii.gz -mask $subj\_mask_upsampled.nii.gz -lmax $lmax +done + +# Generation of 5TT image based on acpc registerred image +5ttgen fsl anatomy/$subj\_t1_acpc.nii.gz $subj\_5tt.nii.gz + +# Generate tractogram + +lmax="8 10 12" +FAs="0.05" +Curv="30 45 60" + +for i in $lmax; do + for j in $FAs; do + for k in $Curv; do + tckgen $subj\_fod_lmax_$lmax\.nii.gz lw_results_r_r_$i\_$j\_$k.tck -number 300 -cutoff $j -angle $k -act lw_5tt.nii.gz -unidirectional -seed_image $subj\_start_right.mif -include $subj\_end_right.mif -stop -step 0.2 -rk4; + tckgen $subj\_fod_lmax_$lmax\.nii.gz lw_results_r_l_$i\_$j\_$k.tck -number 300 -cutoff $j -angle $k -act lw_5tt.nii.gz -unidirectional -seed_image $subj\_start_right.mif -include $subj\_end_left.mif -stop -step 0.2 -rk4; + tckgen $subj\_fod_lmax_$lmax\.nii.gz lw_results_l_r_$i\_$j\_$k.tck -number 300 -cutoff $j -angle $k -act lw_5tt.nii.gz -unidirectional -seed_image $subj\_start_left.mif -include $subj\_end_right.mif -stop -step 0.2 -rk4; + tckgen $subj\_fod_lmax_$lmax\.nii.gz lw_results_l_l_$i\_$j\_$k.tck -number 300 -cutoff $j -angle $k -act lw_5tt.nii.gz -unidirectional -seed_image $subj\_start_left.mif -include $subj\_end_left.mif -stop -step 0.2 -rk4; + echo " " + done + done +done + +altogether="" + +for z in *.tck; do + altogether+="${z} " +done + +tckedit $altogether $subj\_et.tck -force + + +source 03_life.m $subj + +#for r in $subjects; do +# analizuj "$r" +# done + +matlab -nodisplay -r 03_life.m + + diff --git a/02_preparation_for_acpc.sh b/02_preparation_for_acpc.sh new file mode 100644 index 0000000..ca1d232 --- /dev/null +++ b/02_preparation_for_acpc.sh @@ -0,0 +1,21 @@ +ACPC_preparation () { + + # Registration to ACPC space requires preprocessed DWI image in .nii (or .nii.gz) format with corresponding bvecs and bvals. T1 image in ACPC space will be prepared separately + + # Defining data set and its' path + a=${1::-6} + data=/home/auguser2016/dMRI_DATA + high_res=PREPROCESSED_DATA/0.75mm_iso/$a + + # Obtain bvecs and bvals + mrinfo $data/$high_res/$a\_clean_075mm.mif -export_grad_fsl $data/$high_res/$a\_bvecs.bvecs $data/$high_res/$a\_bvals.bvals + + # Transform DW in .mif into .nii.gz + mrconvert $data/$high_res/$a\_clean_075mm.mif $data/$high_res/$a\_clean_075mm.nii.gz +} + +# List of all data_sets used in studies +list=(fe21_1325/ hw91_0844/ ib57_0731/ kw99_0633/ lw37_0977/ nb30_1185/ ow93_0974/ ps94_1516/ rx88_1234/ sj22_1218/ ta14_1065/ tq63_1214/ xs62_1217/ xn78_1085/ uh47_1309/ uf97_1072/ ow93_0974/ ) + +# As during processing of each data set full resources are used, lining them up is an effective approach +for i in ${list[@]}; do ACPC_preparation $i; done diff --git a/02_t1_and_dwi_to_acpc.m b/02_t1_and_dwi_to_acpc.m new file mode 100644 index 0000000..07f4074 --- /dev/null +++ b/02_t1_and_dwi_to_acpc.m @@ -0,0 +1,72 @@ +% Add required folders to path +addpath(genpath('/home/auguser2016/Software/spm12')) +addpath(genpath('/home/auguser2016/Software/vistasoft-master')) +addpath(genpath('/home/auguser2016/Software/mba')) +addpath(genpath('/home/auguser2016/Software/encode')) + +% ACPC registration for T1 in each data set performed manually +% Done manually + +% Perform Diffusion Imaging Processing for all data sets +cd /home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/ +subjects=dir + +for i=3:size(subjects,1) + + id=subjects(i).name + %id='hw91' + % Prepare data in form acceptable by dtiInit + + % ACPC registerred T1 image + anatomy=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/',id,'_t1_acpc.nii.gz') + + % Diffusion Data (motion, eddy and geometrical distortions corrected) + dwi=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_clean_150mm.nii.gz') + + % Bvecs and bvals obtained from corrected DW images + bvecs=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_bvecs.bvecs') + bvals=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_bvals.bvals') + + %bvecs_rot = dlmread(bvecs); + %bvecs_rot(1,:) = -bvecs_rot(1,:); + %dlmwrite(strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id1,'/',id,'_bvecs_rot.bvecs'),bvecs_rot) + %bvecs=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id1,'/',id,'_bvecs_rot.bvecs') + + %bvecs_rot = dlmread(bvecs); + %bvecs_rot(2,:) z= -bvecs_rot(2,:); + %dlmwrite(strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id1,'/',id,'_bvecs_rot.bvecs'),bvecs_rot) + %bvecs=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id1,'/',id,'_bvecs_rot.bvecs') + +% %przypomnienie, że bveci to ważna sprawa i możliwe, że mam filpa na osi y +% bvecs_rot = dlmread(bvecs); +% bvecs_rot(2,:)=-bvecs_rot(2,:); +% dlmwrite(strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_bvecs_rot.bvecs'), bvecs_rot) + + % Output directory + output=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/') + + % Set up dtiInit paramters specific for our study + dwParams = dtiInitParams; + dwParams.rotateBvecsWithRx = 1; + dwParams.rotateBvecsWithCanXform = 1; + dwParams.bvecsFile = fullfile(bvecs); + dwParams.bvalsFile = fullfile(bvals); + dwParams.eddyCorrect = -1; + dwParams.dwOutMm = [0.75 0.75 0.75]; % change this to 1.5 and process 0.75 anyway + %dwParams.dwOutMm = [1.5 1.5 1.5]; + dwParams.phaseEncodeDir = 2; + dwParams.outDir=output + %dwParams.outDir = fullfile(pwd,'dtiInitPreprocessed'); + + % Start Diffusion Imaging preprocessing. + [dt6FileName, outBaseDir] = dtiInit(dwi,... + anatomy, ... + dwParams); + + mrtrix_bfileFromBvecs(strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/',id,'_clean_150mm_aligned_trilin_noMEC.bvecs'),strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/',id,'_clean_150mm_aligned_trilin_noMEC.bvals'),strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/',id,'_mrtrix_grad_table.b')) + + end + + %pathToScript = fullfile('/media/auguser2016/Volume/Test/chiasm','04_postprocess.sh'); + %cmdStr = sprintf('%s', pathToScript); + %system(cmdStr); diff --git a/02_t1_and_dwi_to_acpc.m~ b/02_t1_and_dwi_to_acpc.m~ new file mode 100644 index 0000000..401e4a3 --- /dev/null +++ b/02_t1_and_dwi_to_acpc.m~ @@ -0,0 +1,47 @@ +% Add required folders to path +addpath(genpath('/home/auguser2016/Software/spm12')) +addpath(genpath('/home/auguser2016/Software/vistasoft-master')) +addpath(genpath('/home/auguser2016/Software/mba')) +addpath(genpath('/home/auguser2016/Software/francopestilli-life-aea5433')) + +% ACPC registration for T1 in each data set performed manually +% Done manually + +% Perform Diffusion Imaging Processing for all data sets +cd /home/auguser2016/Projects/Chiasm_RJP_FP/ +subjects=dir + +%for i=3:size(subjects,1) + + %id=subjects(i).name + id='ib57' + + % Prepare data in form acceptable by dtiInit + % ACPC registerred T1 image + anatomy=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP/',id,'/',id,'_t1_acpc.nii.gz') + % Diffusion Data (motion, eddy and geometrical distortions corrected) + dwi=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/',id,'/',id,'_clean_075mm.nii.gz') + % bvecs and bvals obtained from corrected DW images + bvecs=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/',id,'/',id,'.bvecs') + bvals=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/',id,'/',id,'.bvals') + + % Set up dtiInit paramters specific for our study + dwParams = dtiInitParams; + dwParams.rotateBvecsWithRx = 1; + dwParams.rotateBvecsWithCanXform = 1; + dwParams.bvecsFile = fullfile(bvecs); + dwParams.bvalsFile = fullfile(bvals); + dwParams.eddyCorrect = -1; + dwParams.dwOutMm = [1.5 1.5 1.5]; + dwParams.phaseEncodeDir = 2; + dwParams.outDir = fullfile(pwd,'dtiInitPreprocessed'); + + % Start Diffusion Imaging preprocessing. + [dt6FileName, outBaseDir] = dtiInit(dwi,... + anatomy, ... + dwParams); +%end + + + + diff --git a/03_check_and_segment_acpc.sh b/03_check_and_segment_acpc.sh new file mode 100644 index 0000000..c7e7aab --- /dev/null +++ b/03_check_and_segment_acpc.sh @@ -0,0 +1,37 @@ +Check_acpc () { + +subj=$1 + +analysis=/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/$subj + +# Compare coregisterred DW image with anatomy in acpc space +#mrview $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.nii.gz $analysis/$subj\_t1_acpc.nii.gz + +# Compare DW images before and after registration +#mrview $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.nii.gz /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj/$subj\_clean_075mm.mif + +# Compare their bvecs +#gedit $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.bvecs /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj/$subj\_bvecs.bvecs + +# Compare their bvals +#gedit $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.bvals /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj/$subj\_bvals.bvals + +# for ACPC registration +# Perform segmentation of T1 image in ACPC space into WM, GM, subGM, CSF and potential pathological tissue +5ttgen fsl $analysis/$subj\_t1_acpc.nii.gz $analysis/$subj\_t1_acpc_5tt.nii.gz + +# Correct WM segmentation in chiasm - modify opened ROI and close mrview. Modified WM segmentation image will be shown - repeat previous steps until satisfactory result is obtained. Terminate console to finish process +#mrview $analysis/$subj\_t1_acpc_5tt.nii.gz -colourmap 3 -overlay.load $analysis/$subj\_t1_acpc.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $analysis/$subj\_5tt_patch.mif + +#5ttedit $analysis/$subj\_t1_acpc_5tt.nii.gz -wm $analysis/$subj\_5tt_patch.mif $analysis/$subj\_t1_acpc_5tt_modified.nii.gz -force + +#while true; do +# mrview $analysis/$subj\_t1_acpc_5tt_modified.nii.gz -colourmap 3 -overlay.load $analysis/$subj\_t1_acpc.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $analysis/$subj\_5tt_patch.mif +# 5ttedit $analysis/$subj\_t1_acpc_5tt.nii.gz -wm $analysis/$subj\_5tt_patch.mif $analysis/$subj\_t1_acpc_5tt_modified.nii.gz -force +#done + +} + +#for i in /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/*; do Check_acpc ${i: -4}; done + +Check_acpc $1 diff --git a/03_life.m b/03_life.m new file mode 100644 index 0000000..b18b278 --- /dev/null +++ b/03_life.m @@ -0,0 +1,60 @@ + +% Build the file names for the diffusion data, the anatomical MRI. +dwiFile = fullfile('/media/auguser2016/Volume/Test','lw37_aligned_trilin.nii.gz'); +dwiFileRepeat = fullfile('/media/auguser2016/Volume/Test','lw37_aligned_trilin.nii.gz'); +t1File = fullfile('/media/auguser2016/Volume/Test/anatomy','t1_acpc.nii.gz'); +fgFileName = fullfile('/media/auguser2016/Volume/Test','lw_et.tck') + +% The final connectome and data astructure will be saved with this name: +feFileName = 'test_ensemble_tracking'; + +L = 360; % Discretization parameter +fe = feConnectomeInit(dwiFile,fgFileName,feFileName,[],dwiFileRepeat,t1File,L,[1,0]); + +Niter = 500; +fe = feSet(fe,'fit',feFitModel(feGet(fe,'model'),feGet(fe,'dsigdemeaned'),'bbnnls',Niter,'preconditioner')); + +% Remove non zero weighted fascicles +fg = feGet(fe,'fibers acpc'); +w = feGet(fe,'fiber weights'); +positive_w = w > 0; +fg = fgExtract(fg, positive_w, 'keep'); + +% Use MBA to visualize +colors = {[.75 .25 .1]}; +viewCoords = [90,0]; +proportion_to_show = .05; +threshold_length = 10; +slices = {[18 0 0],[0 -40 0],[0 0 -14]}; % Values in mm from AC + +% Prepare the plot of tract +fg = rmfield(fg,'coordspace'); + +% Pick a percentage of fascicles to display (the PN can be too dense for visualization). +%fibs_indx = RandSample(1:length(fg.fibers),round(length(fg.fibers)*proportion_to_show)); +fg.fibers = fg.fibers(1:2000); + +% Display fascicles and anatomy +fig_h = figure('name','Whole Brain','color','k'); +hold on + +% display anatomy +t1 = niftiRead(t1File); +%h = mbaDisplayBrainSlice(t1, slices{1}); +%h = mbaDisplayBrainSlice(t1, slices{2}); + +% Disdplay fasciles +set(gca,'visible','off','ylim',[-108 69],'xlim',[-75 75],'zlim',[-45 78],'Color','w') +[~, light_h] = mbaDisplayConnectome(fg.fibers,fig_h,colors{1},'single'); +delete(light_h) +view(viewCoords(1),viewCoords(2)) +light_h = camlight('right'); +lighting phong; +%set(fig_h,'Units','normalized', 'Position',[0.5 .2 .4 .8]); +set(gcf,'Color',[1 1 1]) +drawnow + +quit + +% save figure to disk +print(fig_h,'my_figure','-dpng','-r600');%'-dpsc2') diff --git a/03_t1_and_dwi_to_acpc (copy).m b/03_t1_and_dwi_to_acpc (copy).m new file mode 100644 index 0000000..1093782 --- /dev/null +++ b/03_t1_and_dwi_to_acpc (copy).m @@ -0,0 +1,54 @@ +% Add required folders to path +addpath(genpath('/home/auguser2016/Software/spm12')) +addpath(genpath('/home/auguser2016/Software/vistasoft-master')) +addpath(genpath('/home/auguser2016/Software/mba')) +addpath(genpath('/home/auguser2016/Software/encode')) + +% ACPC registration for T1 in each data set performed manually +% Done manually + +% Perform Diffusion Imaging Processing for all data sets +cd /home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/ +subjects=dir + +%for i=22:size(subjects,1) + + %id=subjects(i).name + id='uf97' + % Prepare data in form acceptable by dtiInit + + % ACPC registerred T1 image + anatomy=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/',id,'_t1_acpc.nii.gz') + + % Diffusion Data (motion, eddy and geometrical distortions corrected) + dwi=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_clean_150mm.nii.gz') + + % Bvecs and bvals obtained from corrected DW images + bvecs=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_bvecs.bvecs') + bvals=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_bvals.bvals') + + bvecs_rot = dlmread(bvecs); + bvecs_rot(2,:)=-bvecs_rot(2,:); + dlmwrite(strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_bvecs_rot.bvecs'), bvecs_rot) + + % Output directory + output=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'_0p5/') + + % Set up dtiInit paramters specific for our study + dwParams = dtiInitParams; + dwParams.rotateBvecsWithRx = 1; + dwParams.rotateBvecsWithCanXform = 0; + dwParams.bvecsFile = fullfile(bvecs); + dwParams.bvalsFile = fullfile(bvals); + dwParams.eddyCorrect = -1; + dwParams.dwOutMm = [0.75 0.75 0.75]; % change this to 1.5 and process 0.75 anyway + %dwParams.dwOutMm = [1.5 1.5 1.5]; + dwParams.phaseEncodeDir = 2; + dwParams.outDir=output + %dwParams.outDir = fullfile(pwd,'dtiInitPreprocessed'); + + % Start Diffusion Imaging preprocessing. + [dt6FileName, outBaseDir] = dtiInit(dwi,... + anatomy, ... + dwParams); +%end diff --git a/03_t1_and_dwi_to_acpc.m b/03_t1_and_dwi_to_acpc.m new file mode 100644 index 0000000..b1db325 --- /dev/null +++ b/03_t1_and_dwi_to_acpc.m @@ -0,0 +1,51 @@ +% Add required folders to path +addpath(genpath('/home/auguser2016/Software/spm12')) +addpath(genpath('/home/auguser2016/Software/vistasoft-master')) +addpath(genpath('/home/auguser2016/Software/mba')) +addpath(genpath('/home/auguser2016/Software/encode')) + +% ACPC registration for T1 in each data set performed manually +% Done manually + +% Perform Diffusion Imaging Processing for all data sets +cd /home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/ +subjects=dir + +for i=5:size(subjects,1) + + id=subjects(i).name + %id='fe21' + % Prepare data in form acceptable by dtiInit + + % ACPC registerred T1 image + anatomy=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/',id,'_t1_acpc.nii.gz') + + % Diffusion Data (motion, eddy and geometrical distortions corrected) + dwi=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/',id,'/',id,'_clean_075mm.nii.gz') + %dwi=strcat('/media/auguser2016/Volume/Test/lw37.nii.gz') + + % Bvecs and bvals obtained from corrected DW images + bvecs=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/',id,'/',id,'_bvecs.bvecs') + bvals=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/',id,'/',id,'_bvals.bvals') + + % Output directory + output=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/') + + % Set up dtiInit paramters specific for our study + dwParams = dtiInitParams; + dwParams.rotateBvecsWithRx = 1; + dwParams.rotateBvecsWithCanXform = 0; + dwParams.bvecsFile = fullfile(bvecs); + dwParams.bvalsFile = fullfile(bvals); + dwParams.eddyCorrect = -1; + dwParams.dwOutMm = [0.75 0.75 0.75]; % change this to 1.5 and process 0.75 anyway + %dwParams.dwOutMm = [1.5 1.5 1.5]; + dwParams.phaseEncodeDir = 2; + dwParams.outDir=output + %dwParams.outDir = fullfile(pwd,'dtiInitPreprocessed'); + + % Start Diffusion Imaging preprocessing. + [dt6FileName, outBaseDir] = dtiInit(dwi,... + anatomy, ... + dwParams); +end diff --git a/03_working_on_original_lw.m b/03_working_on_original_lw.m new file mode 100644 index 0000000..6ae4082 --- /dev/null +++ b/03_working_on_original_lw.m @@ -0,0 +1,49 @@ +% Add required folders to path +addpath(genpath('/home/auguser2016/Software/spm12')) +addpath(genpath('/home/auguser2016/Software/vistasoft-master')) +addpath(genpath('/home/auguser2016/Software/mba')) +addpath(genpath('/home/auguser2016/Software/encode')) + +% ACPC registration for T1 in each data set performed manually +% Done manually + +% Perform Diffusion Imaging Processing for all data sets +cd /home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/ +subjects=dir + +for i=5:size(subjects,1) + + id=subjects(i).name + % Prepare data in form acceptable by dtiInit + + % ACPC registerred T1 image + anatomy=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/',id,'_t1_acpc.nii.gz') + + % Diffusion Data (motion, eddy and geometrical distortions corrected) + dwi=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_clean_150mm.nii.gz') + + % Bvecs and bvals obtained from corrected DW images + bvecs=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_bvecs.bvecs') + bvals=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_bvals.bvals') + + % Output directory + output=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/') + + % Set up dtiInit paramters specific for our study + dwParams = dtiInitParams; + dwParams.rotateBvecsWithRx = 1; + dwParams.rotateBvecsWithCanXform = 0; + dwParams.bvecsFile = fullfile(bvecs); + dwParams.bvalsFile = fullfile(bvals); + dwParams.eddyCorrect = -1; + dwParams.dwOutMm = [0.75 0.75 0.75]; % change this to 1.5 and process 0.75 anyway + %dwParams.dwOutMm = [1.5 1.5 1.5]; + dwParams.phaseEncodeDir = 2; + dwParams.outDir=output + %dwParams.outDir = fullfile(pwd,'dtiInitPreprocessed'); + + % Start Diffusion Imaging preprocessing. + [dt6FileName, outBaseDir] = dtiInit(dwi,... + anatomy, ... + dwParams); +end diff --git a/04_check_and_segment_acpc.sh b/04_check_and_segment_acpc.sh new file mode 100644 index 0000000..c7e7aab --- /dev/null +++ b/04_check_and_segment_acpc.sh @@ -0,0 +1,37 @@ +Check_acpc () { + +subj=$1 + +analysis=/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/$subj + +# Compare coregisterred DW image with anatomy in acpc space +#mrview $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.nii.gz $analysis/$subj\_t1_acpc.nii.gz + +# Compare DW images before and after registration +#mrview $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.nii.gz /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj/$subj\_clean_075mm.mif + +# Compare their bvecs +#gedit $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.bvecs /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj/$subj\_bvecs.bvecs + +# Compare their bvals +#gedit $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.bvals /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj/$subj\_bvals.bvals + +# for ACPC registration +# Perform segmentation of T1 image in ACPC space into WM, GM, subGM, CSF and potential pathological tissue +5ttgen fsl $analysis/$subj\_t1_acpc.nii.gz $analysis/$subj\_t1_acpc_5tt.nii.gz + +# Correct WM segmentation in chiasm - modify opened ROI and close mrview. Modified WM segmentation image will be shown - repeat previous steps until satisfactory result is obtained. Terminate console to finish process +#mrview $analysis/$subj\_t1_acpc_5tt.nii.gz -colourmap 3 -overlay.load $analysis/$subj\_t1_acpc.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $analysis/$subj\_5tt_patch.mif + +#5ttedit $analysis/$subj\_t1_acpc_5tt.nii.gz -wm $analysis/$subj\_5tt_patch.mif $analysis/$subj\_t1_acpc_5tt_modified.nii.gz -force + +#while true; do +# mrview $analysis/$subj\_t1_acpc_5tt_modified.nii.gz -colourmap 3 -overlay.load $analysis/$subj\_t1_acpc.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $analysis/$subj\_5tt_patch.mif +# 5ttedit $analysis/$subj\_t1_acpc_5tt.nii.gz -wm $analysis/$subj\_5tt_patch.mif $analysis/$subj\_t1_acpc_5tt_modified.nii.gz -force +#done + +} + +#for i in /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/*; do Check_acpc ${i: -4}; done + +Check_acpc $1 diff --git a/04_check_and_segment_acpc.sh~ b/04_check_and_segment_acpc.sh~ new file mode 100644 index 0000000..f80a62c --- /dev/null +++ b/04_check_and_segment_acpc.sh~ @@ -0,0 +1,37 @@ +Check_acpc () { + +subj=$1 + +analysis=/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/$subj + +# Compare coregisterred DW image with anatomy in acpc space +#mrview $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.nii.gz $analysis/$subj\_t1_acpc.nii.gz + +# Compare DW images before and after registration +#mrview $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.nii.gz /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj/$subj\_clean_075mm.mif + +# Compare their bvecs +#gedit $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.bvecs /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj/$subj\_bvecs.bvecs + +# Compare their bvals +#gedit $analysis/$subj\_clean_075mm_aligned_trilin_noMEC.bvals /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj/$subj\_bvals.bvals + +# for ACPC registration +# Perform segmentation of T1 image in ACPC space into WM, GM, subGM, CSF and potential pathological tissue +5ttgen fsl $analysis/$subj\_t1_acpc.nii.gz $analysis/$subj\_t1_acpc_5tt.nii.gz + +# Correct WM segmentation in chiasm - modify opened ROI and close mrview. Modified WM segmentation image will be shown - repeat previous steps until satisfactory result is obtained. Terminate console to finish process +#mrview $analysis/$subj\_t1_acpc_5tt.nii.gz -colourmap 3 -overlay.load $analysis/$subj\_t1_acpc.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $analysis/$subj\_5tt_patch.mif + +#5ttedit $analysis/$subj\_t1_acpc_5tt.nii.gz -wm $analysis/$subj\_5tt_patch.mif $analysis/$subj\_t1_acpc_5tt_modified.nii.gz -force + +#while true; do +# mrview $analysis/$subj\_t1_acpc_5tt_modified.nii.gz -colourmap 3 -overlay.load $analysis/$subj\_t1_acpc.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $analysis/$subj\_5tt_patch.mif +# 5ttedit $analysis/$subj\_t1_acpc_5tt.nii.gz -wm $analysis/$subj\_5tt_patch.mif $analysis/$subj\_t1_acpc_5tt_modified.nii.gz -force +#done + +} + +#for i in /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/*; do Check_acpc ${i: -4}; done + +Check_acpc hw91 diff --git a/04_postprocess.sh b/04_postprocess.sh new file mode 100644 index 0000000..a1cf419 --- /dev/null +++ b/04_postprocess.sh @@ -0,0 +1,62 @@ + #!/bin/bash + +General_postprocessing () { + +subj=${1} + +analysis=/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/$subj +#analysis=/home/auguser2016/Projects/hw91_canxform + +dwi2mask $analysis/$subj\_clean_150mm_aligned_trilin_noMEC.nii.gz $analysis/$subj\_acpc_mask.nii.gz -grad $analysis/$subj\_mrtrix_grad_table.b + +# Estimating single fibre response (SFR) for lmax=6 +dwi2response tournier $analysis/$subj\_clean_150mm_aligned_trilin_noMEC.nii.gz $analysis/$subj\_acpc_SFR.txt -shell 1600 -lmax 6 -mask $analysis/$subj\_acpc_mask.nii.gz -voxels $analysis/$subj\_acpc_SFR_voxels.nii.gz -grad $analysis/$subj\_mrtrix_grad_table.b + +# Estimation of Fiber Orientation Distribution function for 3 lmax +for lmax in 8 10 12; do + dwiextract $analysis/$subj\_clean_150mm_aligned_trilin_noMEC.nii.gz -grad $analysis/$subj\_mrtrix_grad_table.b - | dwi2fod msmt_csd - $analysis/$subj\_acpc_SFR.txt -lmax $lmax $analysis/$subj\_acpc_FOD_lmax_$lmax.mif -mask $analysis/$subj\_acpc_mask.nii.gz +done +} + +current=$(pwd) + +cd /home/auguser2016/Projects/Chiasm_RJP_FP_ACPC +for i in */; do + General_postprocessing ${i::-1}; +done + +cd $current + + +''' +# List of all data_sets used in studies +list=(fe21_1325/ hw91_0844/ ib57_0731/ kw99_0633/ lw37_0977/ nb30_1185/ ow93_0974/ ps94_1516/ rx88_1234/ sj22_1218/ ta14_1065/ tq63_1214/ xs62_1217/ xn78_1085/ uh47_1309/ uf97_1072/ ow93_0974/ ) + +list=(hw91_some/) #ow93_0974/ ) + + +# This way 4 sets are processed at once - we can't process all at once due to RAM usage, however preprocessing one at a time is ineffective because e.g. eddy from dwipreproc is using only one core. There is a version of eddy operating on multicore processors, but it's not running on this machine + +for i in ${list[@]:0:3}; do General_postprocessing $i; done +wait +for i in ${list[@]:3:6}; do General_postprocessing $i & done +wait +for i in ${list[@]:6:9}; do General_postprocessing $i & done +wait +for i in ${list[@]:9:12}; do General_postprocessing $i & done +wait +for i in ${list[@]:12:15}; do General_postprocessing $i & done +wait +for i in ${list[@]:15:18}; do General_postprocessing $i & done +wait +for i in ${list[@]:18:21}; do General_postprocessing $i & done +wait +for i in ${list[@]:21:24}; do General_postprocessing $i & done +wait +for i in ${list[@]:24:27}; do General_postprocessing $i & done +wait +for i in ${list[@]:27:30}; do General_postprocessing $i & done +wait +for i in ${list[@]:30:33}; do General_postprocessing $i & done +wait +''' diff --git a/04_postprocess.sh~ b/04_postprocess.sh~ new file mode 100644 index 0000000..a1cf419 --- /dev/null +++ b/04_postprocess.sh~ @@ -0,0 +1,62 @@ + #!/bin/bash + +General_postprocessing () { + +subj=${1} + +analysis=/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/$subj +#analysis=/home/auguser2016/Projects/hw91_canxform + +dwi2mask $analysis/$subj\_clean_150mm_aligned_trilin_noMEC.nii.gz $analysis/$subj\_acpc_mask.nii.gz -grad $analysis/$subj\_mrtrix_grad_table.b + +# Estimating single fibre response (SFR) for lmax=6 +dwi2response tournier $analysis/$subj\_clean_150mm_aligned_trilin_noMEC.nii.gz $analysis/$subj\_acpc_SFR.txt -shell 1600 -lmax 6 -mask $analysis/$subj\_acpc_mask.nii.gz -voxels $analysis/$subj\_acpc_SFR_voxels.nii.gz -grad $analysis/$subj\_mrtrix_grad_table.b + +# Estimation of Fiber Orientation Distribution function for 3 lmax +for lmax in 8 10 12; do + dwiextract $analysis/$subj\_clean_150mm_aligned_trilin_noMEC.nii.gz -grad $analysis/$subj\_mrtrix_grad_table.b - | dwi2fod msmt_csd - $analysis/$subj\_acpc_SFR.txt -lmax $lmax $analysis/$subj\_acpc_FOD_lmax_$lmax.mif -mask $analysis/$subj\_acpc_mask.nii.gz +done +} + +current=$(pwd) + +cd /home/auguser2016/Projects/Chiasm_RJP_FP_ACPC +for i in */; do + General_postprocessing ${i::-1}; +done + +cd $current + + +''' +# List of all data_sets used in studies +list=(fe21_1325/ hw91_0844/ ib57_0731/ kw99_0633/ lw37_0977/ nb30_1185/ ow93_0974/ ps94_1516/ rx88_1234/ sj22_1218/ ta14_1065/ tq63_1214/ xs62_1217/ xn78_1085/ uh47_1309/ uf97_1072/ ow93_0974/ ) + +list=(hw91_some/) #ow93_0974/ ) + + +# This way 4 sets are processed at once - we can't process all at once due to RAM usage, however preprocessing one at a time is ineffective because e.g. eddy from dwipreproc is using only one core. There is a version of eddy operating on multicore processors, but it's not running on this machine + +for i in ${list[@]:0:3}; do General_postprocessing $i; done +wait +for i in ${list[@]:3:6}; do General_postprocessing $i & done +wait +for i in ${list[@]:6:9}; do General_postprocessing $i & done +wait +for i in ${list[@]:9:12}; do General_postprocessing $i & done +wait +for i in ${list[@]:12:15}; do General_postprocessing $i & done +wait +for i in ${list[@]:15:18}; do General_postprocessing $i & done +wait +for i in ${list[@]:18:21}; do General_postprocessing $i & done +wait +for i in ${list[@]:21:24}; do General_postprocessing $i & done +wait +for i in ${list[@]:24:27}; do General_postprocessing $i & done +wait +for i in ${list[@]:27:30}; do General_postprocessing $i & done +wait +for i in ${list[@]:30:33}; do General_postprocessing $i & done +wait +''' diff --git a/04_psychophysics.sh~ b/04_psychophysics.sh~ new file mode 100644 index 0000000..e69de29 diff --git a/05_ensemble_tracking.txt b/05_ensemble_tracking.txt new file mode 100644 index 0000000..754fb82 --- /dev/null +++ b/05_ensemble_tracking.txt @@ -0,0 +1,258 @@ +# 1 Argument: ID of participant + +# Requirements: +# Data preprocessed +# DWI corregisterd to anatomy +# Done and corrected 5TT segmentation + +ROIs_tracking_anatomy () { + +subj=$1 +data=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj +output=/home/auguser2016/Projects/Chiasm_RJP_FP/$subj +mkdir $output +mkdir $output/Tractography_dwi +cd $output + +# Manually create 4 ROIs for given data set +#mrview $data/$subj\_clean_075mm_FOD.mif -overlay.load $data/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz + +# SFR Estimation +lmax=6 #This should be as low as possible, optimal for most participants value #SFR for lower lmax, keep it and use higher 8,10,12 for fod estimation OR Get SFR with optimal lmax for most participants and use given lmax for all participants. Lower lmax preffered +#dwi2response tournier $data/$subj\_clean_075mm.mif -shell 1600 -lmax $lmax -mask $data/$subj\_clean_075mm_mask.mif $output/$subj\_sfr.txt -voxels $data/Temporal/$subj\_SFR_voxels.mif + +# Estimation of Fiber Orientation Distribution function +#for lmax in 8 10 12; do +# dwiextract $data/$subj\_clean_075mm.mif - | dwi2fod msmt_csd - $output/$subj\_sfr.txt -lmax $lmax $output/$subj\_FOD_lmax_$lmax.mif -mask $data/$subj\_clean_075mm_mask.mif +#done + +# Generation of 5TT image based on acpc registerred image +#5ttgen fsl anatomy/$subj\_t1_acpc.nii.gz $subj\_5tt.nii.gz + +# Generate tractogram + +lmax="8 10 12" +FAs="0.05" +Curv="30 45 60" + +tt_image=$data/$subj\_5tt_for_act_modified.nii.gz + +for i in $lmax; do + for j in $FAs; do + for k in $Curv; do + echo $subj\_results_r_r_$i\_$j\_$k + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_r_l_$i\_$j\_$k + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_l_r_$i\_$j\_$k + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_l_l_$i\_$j\_$k + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 100000; + + #echo $subj\_results_r_r_$i\_$j\_$k\_no_ACT + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_r_l_$i\_$j\_$k\_no_ACT + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_l_r_$i\_$j\_$k\_no_ACT + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_l_l_$i\_$j\_$k\_no_ACT + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sl.mif -include $output//$subj\_el.mif -stop -step 0.2 -maxnum 100000; + + + echo $subj\_results_r_r_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $data/$subj\_5tt_for_act_modified.nii.gz -unidirectional -seed_image $output/$subj\_start_right.mif -include $output/$subj\_end_right.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_r_l_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $data/$subj\_5tt_for_act_modified.nii.gz -unidirectional -seed_image $output/$subj\_start_right.mif -include $output/$subj\_end_left.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_l_r_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $data/$subj\_5tt_for_act_modified.nii.gz -unidirectional -seed_image $output/$subj\_start_left.mif -include $output/$subj\_end_right.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_l_l_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $data/$subj\_5tt_for_act_modified.nii.gz -unidirectional -seed_image $output/$subj\_start_left.mif -include $output/$subj\_end_left.mif -stop -step 0.2 -maxnum 100000; + + echo $subj\_results_r_r_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff 0.10 -angle $k -unidirectional -seed_image $output/$subj\_start_right.mif -include $output/$subj\_end_right.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_r_l_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff 0.10 -angle $k -unidirectional -seed_image $output/$subj\_start_right.mif -include $output/$subj\_end_left.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_l_r_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff 0.10 -angle $k -unidirectional -seed_image $output/$subj\_start_left.mif -include $output/$subj\_end_right.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_l_l_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff 0.10 -angle $k -unidirectional -seed_image $output/$subj\_start_left.mif -include $output//$subj\_end_left.mif -stop -step 0.2 -maxnum 100000; + + echo " " + done + done +done + + + +cd $output/Tractography + +tckedit *_??.tck $subj\_ACT.tck +tckedit *r_r*_??.tck $subj\_r_r_ACT.tck +tckedit *r_l*_??.tck $subj\_r_l_ACT.tck +tckedit *l_r*_??.tck $subj\_l_r_ACT.tck +tckedit *l_l*_??.tck $subj\_l_l_ACT.tck + +tckedit *_no_ACT.tck $subj\_no_ACT.tck +tckedit *r_r*_no_ACT.tck $subj\_r_r_no_ACT.tck +tckedit *r_l*_no_ACT.tck $subj\_r_l_no_ACT.tck +tckedit *l_r*_no_ACT.tck $subj\_l_r_no_ACT.tck +tckedit *l_l*_no_ACT.tck $subj\_l_l_no_ACT.tck + +} + +ROIs_tracking_dwi () { + +subj=$1 +data=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj +output=/home/auguser2016/Projects/Chiasm_RJP_FP/$subj +mkdir $output +mkdir $output/Tractography_dwi +cd $output + +# Manually create 4 ROIs for given data set +#mrview $data/$subj\_clean_075mm_FOD.mif -overlay.load $data/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz + +# SFR Estimation +lmax=6 #This should be as low as possible, optimal for most participants value #SFR for lower lmax, keep it and use higher 8,10,12 for fod estimation OR Get SFR with optimal lmax for most participants and use given lmax for all participants. Lower lmax preffered +#dwi2response tournier $data/$subj\_clean_075mm.mif -shell 1600 -lmax $lmax -mask $data/$subj\_clean_075mm_mask.mif $output/$subj\_sfr.txt -voxels $data/Temporal/$subj\_SFR_voxels.mif + +# Estimation of Fiber Orientation Distribution function +#for lmax in 8 10 12; do +# dwiextract $data/$subj\_clean_075mm.mif - | dwi2fod msmt_csd - $output/$subj\_sfr.txt -lmax $lmax $output/$subj\_FOD_lmax_$lmax.mif -mask $data/$subj\_clean_075mm_mask.mif +#done + +# Generation of 5TT image based on acpc registerred image +#5ttgen fsl anatomy/$subj\_t1_acpc.nii.gz $subj\_5tt.nii.gz + +# Generate tractogram + +lmax="8 10 12" +FAs="0.05" +Curv="30 45 60" + +tt_image=$data/$subj\_5tt_for_act_modified_dwi.nii.gz + +for i in $lmax; do + for j in $FAs; do + for k in $Curv; do + echo $subj\_results_r_r_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_r_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_r_l_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_r_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_l_r_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_l_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_l_l_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_l_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 10000000; + + echo $subj\_results_r_r_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_r_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_r_l_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_r_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_l_r_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_l_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_l_l_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_l_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sl.mif -include $output//$subj\_el.mif -stop -step 0.2 -maxnum 10000000; + + echo " " + done + done +done + + + +cd $output/Tractography_dwi + +tckedit *_??.tck $subj\_ACT.tck +tckedit *r_r*_??.tck $subj\_r_r_ACT.tck +tckedit *r_l*_??.tck $subj\_r_l_ACT.tck +tckedit *l_r*_??.tck $subj\_l_r_ACT.tck +tckedit *l_l*_??.tck $subj\_l_l_ACT.tck + +tckedit *_no_ACT.tck $subj\_no_ACT.tck +tckedit *r_r*_no_ACT.tck $subj\_r_r_no_ACT.tck +tckedit *r_l*_no_ACT.tck $subj\_r_l_no_ACT.tck +tckedit *l_r*_no_ACT.tck $subj\_l_r_no_ACT.tck +tckedit *l_l*_no_ACT.tck $subj\_l_l_no_ACT.tck + +} + + +ACPC_tracking () { + +subj=$1 +data=/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/$subj +output=/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/$subj + + +mkdir $output/Tractography +cd $output + +# SFR Estimation +#lmax=6 #This should be as low as possible, optimal for most participants value #SFR for lower lmax, keep it and use higher 8,10,12 for fod estimation OR Get SFR with optimal lmax for most participants and use given lmax for all participants. Lower lmax preffered +#dwi2response tournier $data/$subj\_clean_075mm.mif -shell 1600 -lmax $lmax -mask $data/$subj\_clean_075mm_mask.mif $output/$subj\_sfr.txt -voxels $data/Temporal/$subj\_SFR_voxels.mif + +# Estimation of Fiber Orientation Distribution function +#for lmax in 8 10 12; do +# dwiextract $data/$subj\_clean_075mm.mif - | dwi2fod msmt_csd - $output/$subj\_sfr.txt -lmax $lmax $output/$subj\_FOD_lmax_$lmax.mif -mask $data/$subj\_clean_075mm_mask.mif +#done + +# Generation of 5TT image based on acpc registerred image +#5ttgen fsl anatomy/$subj\_t1_acpc.nii.gz $subj\_5tt.nii.gz + +# Generate tractogram + +lmax="8 10 12" +FAs="0.05" +Curv="30 45 60" + +tt_image=$data/$subj\_t1_acpc_5tt_corrected.nii.gz +#tt_image=$data/$subj\_t1_acpc_5tt.nii.gz + +for i in $lmax; do + for j in $FAs; do + for k in $Curv; do + echo $subj\_results_r_r_$i\_$j\_$k + tckgen $output/$subj\_acpc_FOD_lmax_$i.mif Tractography/$subj\_results_r_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/ROIs/$subj\_sr.mif -include $output/ROIs/$subj\_er.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_r_l_$i\_$j\_$k + tckgen $output/$subj\_acpc_FOD_lmax_$i.mif Tractography/$subj\_results_r_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/ROIs/$subj\_sr.mif -include $output/ROIs/$subj\_el.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_l_r_$i\_$j\_$k + tckgen $output/$subj\_acpc_FOD_lmax_$i.mif Tractography/$subj\_results_l_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/ROIs/$subj\_sl.mif -include $output/ROIs/$subj\_er.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_l_l_$i\_$j\_$k + tckgen $output/$subj\_acpc_FOD_lmax_$i.mif Tractography/$subj\_results_l_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/ROIs/$subj\_sl.mif -include $output/ROIs/$subj\_el.mif -stop -step 0.2 -maxnum 10000000; + + echo " " + done + done +done + + + +cd $output/Tractography + +tckedit *_??.tck $subj\_ACT.tck +tckedit *r_r*_??.tck $subj\_r_r_ACT.tck +tckedit *r_l*_??.tck $subj\_r_l_ACT.tck +tckedit *l_r*_??.tck $subj\_l_r_ACT.tck +tckedit *l_l*_??.tck $subj\_l_l_ACT.tck + +} + + + + + +#for i in /home/auguser2016/Projects/Chiasm_RJP_FP/*; do ROIs_tracking ${i: -4}; done + + + +ACPC_tracking $1 + +#source 03_life.m $subj + +#for r in $subjects; do +# analizuj "$r" +# done + +#matlab -nodisplay -r 03_life.m + + diff --git a/05_ensemble_tracking.txt~ b/05_ensemble_tracking.txt~ new file mode 100644 index 0000000..26238c5 --- /dev/null +++ b/05_ensemble_tracking.txt~ @@ -0,0 +1,198 @@ +# 1 Argument: ID of participant + +# Requirements: +# Data preprocessed +# DWI corregisterd to anatomy +# Done and corrected 5TT segmentation + +ROIs_tracking_anatomy () { + +subj=$1 +data=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj +output=/home/auguser2016/Projects/Chiasm_RJP_FP/$subj +mkdir $output +mkdir $output/Tractography_dwi +cd $output + +# Manually create 4 ROIs for given data set +#mrview $data/$subj\_clean_075mm_FOD.mif -overlay.load $data/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz + +# SFR Estimation +lmax=6 #This should be as low as possible, optimal for most participants value #SFR for lower lmax, keep it and use higher 8,10,12 for fod estimation OR Get SFR with optimal lmax for most participants and use given lmax for all participants. Lower lmax preffered +#dwi2response tournier $data/$subj\_clean_075mm.mif -shell 1600 -lmax $lmax -mask $data/$subj\_clean_075mm_mask.mif $output/$subj\_sfr.txt -voxels $data/Temporal/$subj\_SFR_voxels.mif + +# Estimation of Fiber Orientation Distribution function +#for lmax in 8 10 12; do +# dwiextract $data/$subj\_clean_075mm.mif - | dwi2fod msmt_csd - $output/$subj\_sfr.txt -lmax $lmax $output/$subj\_FOD_lmax_$lmax.mif -mask $data/$subj\_clean_075mm_mask.mif +#done + +# Generation of 5TT image based on acpc registerred image +#5ttgen fsl anatomy/$subj\_t1_acpc.nii.gz $subj\_5tt.nii.gz + +# Generate tractogram + +lmax="8 10 12" +FAs="0.05" +Curv="30 45 60" + +tt_image=$data/$subj\_5tt_for_act_modified.nii.gz + +for i in $lmax; do + for j in $FAs; do + for k in $Curv; do + echo $subj\_results_r_r_$i\_$j\_$k + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_r_l_$i\_$j\_$k + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_l_r_$i\_$j\_$k + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_l_l_$i\_$j\_$k + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 100000; + + #echo $subj\_results_r_r_$i\_$j\_$k\_no_ACT + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_r_l_$i\_$j\_$k\_no_ACT + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_l_r_$i\_$j\_$k\_no_ACT + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 100000; + #echo $subj\_results_l_l_$i\_$j\_$k\_no_ACT + #tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sl.mif -include $output//$subj\_el.mif -stop -step 0.2 -maxnum 100000; + + + echo $subj\_results_r_r_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $data/$subj\_5tt_for_act_modified.nii.gz -unidirectional -seed_image $output/$subj\_start_right.mif -include $output/$subj\_end_right.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_r_l_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $data/$subj\_5tt_for_act_modified.nii.gz -unidirectional -seed_image $output/$subj\_start_right.mif -include $output/$subj\_end_left.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_l_r_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $data/$subj\_5tt_for_act_modified.nii.gz -unidirectional -seed_image $output/$subj\_start_left.mif -include $output/$subj\_end_right.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_l_l_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $data/$subj\_5tt_for_act_modified.nii.gz -unidirectional -seed_image $output/$subj\_start_left.mif -include $output/$subj\_end_left.mif -stop -step 0.2 -maxnum 100000; + + echo $subj\_results_r_r_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff 0.10 -angle $k -unidirectional -seed_image $output/$subj\_start_right.mif -include $output/$subj\_end_right.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_r_l_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_r_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff 0.10 -angle $k -unidirectional -seed_image $output/$subj\_start_right.mif -include $output/$subj\_end_left.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_l_r_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff 0.10 -angle $k -unidirectional -seed_image $output/$subj\_start_left.mif -include $output/$subj\_end_right.mif -stop -step 0.2 -maxnum 100000; + echo $subj\_results_l_l_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography/$subj\_results_l_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff 0.10 -angle $k -unidirectional -seed_image $output/$subj\_start_left.mif -include $output//$subj\_end_left.mif -stop -step 0.2 -maxnum 100000; + + echo " " + done + done +done + + + +cd $output/Tractography + +tckedit *_??.tck $subj\_ACT.tck +tckedit *r_r*_??.tck $subj\_r_r_ACT.tck +tckedit *r_l*_??.tck $subj\_r_l_ACT.tck +tckedit *l_r*_??.tck $subj\_l_r_ACT.tck +tckedit *l_l*_??.tck $subj\_l_l_ACT.tck + +tckedit *_no_ACT.tck $subj\_no_ACT.tck +tckedit *r_r*_no_ACT.tck $subj\_r_r_no_ACT.tck +tckedit *r_l*_no_ACT.tck $subj\_r_l_no_ACT.tck +tckedit *l_r*_no_ACT.tck $subj\_l_r_no_ACT.tck +tckedit *l_l*_no_ACT.tck $subj\_l_l_no_ACT.tck + +} + +ROIs_tracking_dwi () { + +subj=$1 +data=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj +output=/home/auguser2016/Projects/Chiasm_RJP_FP/$subj +mkdir $output +mkdir $output/Tractography_dwi +cd $output + +# Manually create 4 ROIs for given data set +#mrview $data/$subj\_clean_075mm_FOD.mif -overlay.load $data/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz + +# SFR Estimation +lmax=6 #This should be as low as possible, optimal for most participants value #SFR for lower lmax, keep it and use higher 8,10,12 for fod estimation OR Get SFR with optimal lmax for most participants and use given lmax for all participants. Lower lmax preffered +#dwi2response tournier $data/$subj\_clean_075mm.mif -shell 1600 -lmax $lmax -mask $data/$subj\_clean_075mm_mask.mif $output/$subj\_sfr.txt -voxels $data/Temporal/$subj\_SFR_voxels.mif + +# Estimation of Fiber Orientation Distribution function +#for lmax in 8 10 12; do +# dwiextract $data/$subj\_clean_075mm.mif - | dwi2fod msmt_csd - $output/$subj\_sfr.txt -lmax $lmax $output/$subj\_FOD_lmax_$lmax.mif -mask $data/$subj\_clean_075mm_mask.mif +#done + +# Generation of 5TT image based on acpc registerred image +#5ttgen fsl anatomy/$subj\_t1_acpc.nii.gz $subj\_5tt.nii.gz + +# Generate tractogram + +lmax="8 10 12" +FAs="0.05" +Curv="30 45 60" + +tt_image=$data/$subj\_5tt_for_act_modified_dwi.nii.gz + +for i in $lmax; do + for j in $FAs; do + for k in $Curv; do + echo $subj\_results_r_r_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_r_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_r_l_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_r_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_l_r_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_l_r_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_l_l_$i\_$j\_$k + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_l_l_$i\_$j\_$k.tck -number 200 -cutoff $j -angle $k -act $tt_image -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 10000000; + + echo $subj\_results_r_r_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_r_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_r_l_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_r_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sr.mif -include $output/$subj\_el.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_l_r_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_l_r_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sl.mif -include $output/$subj\_er.mif -stop -step 0.2 -maxnum 10000000; + echo $subj\_results_l_l_$i\_$j\_$k\_no_ACT + tckgen $output/$subj\_FOD_lmax_$i.mif Tractography_dwi/$subj\_results_l_l_$i\_$j\_$k\_no_ACT.tck -number 200 -cutoff $j -angle $k -unidirectional -seed_image $output/$subj\_sl.mif -include $output//$subj\_el.mif -stop -step 0.2 -maxnum 10000000; + + echo " " + done + done +done + + + +cd $output/Tractography_dwi + +tckedit *_??.tck $subj\_ACT.tck +tckedit *r_r*_??.tck $subj\_r_r_ACT.tck +tckedit *r_l*_??.tck $subj\_r_l_ACT.tck +tckedit *l_r*_??.tck $subj\_l_r_ACT.tck +tckedit *l_l*_??.tck $subj\_l_l_ACT.tck + +tckedit *_no_ACT.tck $subj\_no_ACT.tck +tckedit *r_r*_no_ACT.tck $subj\_r_r_no_ACT.tck +tckedit *r_l*_no_ACT.tck $subj\_r_l_no_ACT.tck +tckedit *l_r*_no_ACT.tck $subj\_l_r_no_ACT.tck +tckedit *l_l*_no_ACT.tck $subj\_l_l_no_ACT.tck + +} + +#for i in /home/auguser2016/Projects/Chiasm_RJP_FP/*; do ROIs_tracking ${i: -4}; done + + +ROIs_tracking_dwi fe21 +ROIs_tracking_dwi hw91 +ROIs_tracking_dwi ib57 +ROIs_tracking_dwi kw99 +ROIs_tracking_dwi lw37 +ROIs_tracking_dwi nb30 +ROIs_tracking_dwi ow97 + +#source 03_life.m $subj + +#for r in $subjects; do +# analizuj "$r" +# done + +#matlab -nodisplay -r 03_life.m + + diff --git a/06_life.m b/06_life.m new file mode 100644 index 0000000..5f56da8 --- /dev/null +++ b/06_life.m @@ -0,0 +1,59 @@ +% Build the file names for the diffusion data, the anatomical MRI. +dwiFile = fullfile('/media/auguser2016/Volume/Test','lw37_aligned_trilin.nii.gz'); +dwiFileRepeat = fullfile('/media/auguser2016/Volume/Test','lw37_aligned_trilin.nii.gz'); +t1File = fullfile('/media/auguser2016/Volume/Test/anatomy','t1_acpc.nii.gz'); +fgFileName = fullfile('/media/auguser2016/Volume/Test','lw_et.tck') + +% The final connectome and data astructure will be saved with this name: +feFileName = 'test_ensemble_tracking'; + +L = 360; % Discretization parameter +fe = feConnectomeInit(dwiFile,fgFileName,feFileName,[],dwiFileRepeat,t1File,L,[0,1]); + +Niter = 500; +fe = feSet(fe,'fit',feFitModel(feGet(fe,'model'),feGet(fe,'dsigdemeaned'),'bbnnls',Niter,'preconditioner')); + +% Remove non zero weighted fascicles +fg = feGet(fe,'fibers acpc'); +w = feGet(fe,'fiber weights'); +positive_w = w > 0; +fg = fgExtract(fg, positive_w, 'keep'); + +% Use MBA to visualize +colors = {[.75 .25 .1]}; +viewCoords = [90,0]; +proportion_to_show = .05; +threshold_length = 10; +slices = {[18 0 0],[0 -40 0],[0 0 -14]}; % Values in mm from AC + +% Prepare the plot of tract +fg = rmfield(fg,'coordspace'); + +% Pick a percentage of fascicles to display (the PN can be too dense for visualization). +%fibs_indx = RandSample(1:length(fg.fibers),round(length(fg.fibers)*proportion_to_show)); +fg.fibers = fg.fibers(1:2000); + +% Display fascicles and anatomy +fig_h = figure('name','Whole Brain','color','k'); +hold on + +% display anatomy +t1 = niftiRead(t1File); +%h = mbaDisplayBrainSlice(t1, slices{1}); +%h = mbaDisplayBrainSlice(t1, slices{2}); + +% Disdplay fasciles +set(gca,'visible','off','ylim',[-108 69],'xlim',[-75 75],'zlim',[-45 78],'Color','w') +[~, light_h] = mbaDisplayConnectome(fg.fibers,fig_h,colors{1},'single'); +delete(light_h) +view(viewCoords(1),viewCoords(2)) +light_h = camlight('right'); +lighting phong; +%set(fig_h,'Units','normalized', 'Position',[0.5 .2 .4 .8]); +set(gcf,'Color',[1 1 1]) +drawnow + +quit + +% save figure to disk +print(fig_h,'my_figure','-dpng','-r600');%'-dpsc2') diff --git a/06_life_no_ACPC.m b/06_life_no_ACPC.m new file mode 100644 index 0000000..47dbe01 --- /dev/null +++ b/06_life_no_ACPC.m @@ -0,0 +1,79 @@ +% Add required folders to path +addpath(genpath('/home/auguser2016/Software/spm12')) +addpath(genpath('/home/auguser2016/Software/vistasoft-master')) +addpath(genpath('/home/auguser2016/Software/mba')) +addpath(genpath('/home/auguser2016/Software/encode')) + +% Define participant +subj ='kw99' + +% Define relevant folders +general_data=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/',subj,'/') +specific_data=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP/',subj,'/') +current = pwd + +cd(general_data) + +% Build the file names for the diffusion data, the anatomical MRI. +dwiFile = fullfile(strcat(general_data,subj,'_clean_075mm.nii.gz')); +dwiFileRepeat = fullfile(strcat(general_data,subj,'_clean_075mm.nii.gz')); +t1File = fullfile(strcat(general_data,subj,'_075mm_anatomical_coreg2nodif_undist_2.nii.gz')); +fgFileName = fullfile(strcat(specific_data,'Tractography_dwi/',subj,'_ACT.tck')); +%fgFileName = fullfile(strcat(specific_data,'Tractography/',subj,'_good.tck')); + +% The final connectome and data astructure will be saved with this name: +feFileName = 'test_ensemble_tracking'; + +L = 360; % Discretization parameter +fe = feConnectomeInit(dwiFile,fgFileName,feFileName,[],dwiFileRepeat,t1File,L,[0,1]); + +Niter = 500; +fe = feSet(fe,'fit',feFitModel(feGet(fe,'model'),feGet(fe,'dsigdemeaned'),'bbnnls',Niter,'preconditioner')); + +% Remove non zero weighted fascicles +fg = feGet(fe,'fibers acpc'); +w = feGet(fe,'fiber weights'); +positive_w = w > 0; +fg = fgExtract(fg, positive_w, 'keep'); + +cd(specific_data) +fgWrite(fg,strcat(subj,'_dwi_filtered.mat')) + +% Use MBA to visualize +colors = {[.75 .25 .1]}; +viewCoords = [90,0]; +proportion_to_show = .05; +threshold_length = 10; +slices = {[18 0 0],[0 -40 0],[0 0 -14]}; % Values in mm from AC + +% Prepare the plot of tract +fg = rmfield(fg,'coordspace'); + +% Pick a percentage of fascicles to display (the PN can be too dense for visualization). +%fibs_indx = RandSample(1:length(fg.fibers),round(length(fg.fibers)*proportion_to_show)); +%fg.fibers = fg.fibers(1:1000); + +% Display fascicles and anatomy +fig_h = figure('name','Whole Brain','color','k'); +hold on + +% display anatomy +%t1 = niftiRead(t1File); +%h = mbaDisplayBrainSlice(t1, slices{1}); +%h = mbaDisplayBrainSlice(t1, slices{2}); + +% Display fasciles +set(gca,'visible','off','ylim',[-108 69],'xlim',[-75 75],'zlim',[-45 78],'Color','w') +[~, light_h] = mbaDisplayConnectome(fg.fibers,fig_h,colors{1},'single'); +delete(light_h) +view(viewCoords(1),viewCoords(2)) +light_h = camlight('right'); +lighting phong; +%set(fig_h,'Units','normalized', 'Position',[0.5 .2 .4 .8]); +set(gcf,'Color',[1 1 1]) +drawnow + +%quit + +% save figure to disk +%print(fig_h,'my_figure','-dpng','-r600');%'-dpsc2') \ No newline at end of file diff --git a/07_results_evaluation.txt b/07_results_evaluation.txt new file mode 100644 index 0000000..5cdabdc --- /dev/null +++ b/07_results_evaluation.txt @@ -0,0 +1,59 @@ +Evaluate_anatomy (){ + +cur_path=$(pwd) +subj=$1 + +analysis_folder=/home/auguser2016/Projects/Chiasm_RJP_FP/$subj +mkdir $analysis_folder/ROIs +cd $analysis_folder/ROIs + +mrconvert $analysis_folder/$subj\_start_right.mif $analysis_folder/ROIs/$subj\_start_right.nii +mrconvert $analysis_folder/$subj\_start_left.mif $analysis_folder/ROIs/$subj\_start_left.nii +mrconvert $analysis_folder/$subj\_end_right.mif $analysis_folder/ROIs/$subj\_end_right.nii +mrconvert $analysis_folder/$subj\_end_left.mif $analysis_folder/ROIs/$subj\_end_left.nii + +matlab -r "addpath(genpath('/home/auguser2016/Software/vistasoft-master'));addpath(genpath('/home/auguser2016/Software/mba'));addpath(genpath('/home/auguser2016/Software/encode'));dtiImportRoiFromNifti('$analysis_folder/ROIs/${subj}_start_right.nii','${subj}_start_right');dtiImportRoiFromNifti('$analysis_folder/ROIs/${subj}_start_left.nii','${subj}_start_left');dtiImportRoiFromNifti('$analysis_folder/ROIs/${subj}_end_right.nii','${subj}_end_right');dtiImportRoiFromNifti('$analysis_folder/ROIs/${subj}_end_left.nii','${subj}_end_left')" + +cd $cur_path +} + +Evaluate_dwi (){ + +cur_path=$(pwd) +subj=$1 + +analysis_folder=/home/auguser2016/Projects/Chiasm_RJP_FP/$subj +mkdir $analysis_folder/ROIs_dwi +cd $analysis_folder/ROIs_dwi + +mrconvert $analysis_folder/$subj\_sr.mif $analysis_folder/ROIs_dwi/$subj\_start_right.nii +mrconvert $analysis_folder/$subj\_sl.mif $analysis_folder/ROIs_dwi/$subj\_start_left.nii +mrconvert $analysis_folder/$subj\_er.mif $analysis_folder/ROIs_dwi/$subj\_end_right.nii +mrconvert $analysis_folder/$subj\_el.mif $analysis_folder/ROIs_dwi/$subj\_end_left.nii + +matlab -r "addpath(genpath('/home/auguser2016/Software/vistasoft-master'));addpath(genpath('/home/auguser2016/Software/mba'));addpath(genpath('/home/auguser2016/Software/encode'));dtiImportRoiFromNifti('$analysis_folder/ROIs_dwi/${subj}_start_right.nii','${subj}_start_right');dtiImportRoiFromNifti('$analysis_folder/ROIs_dwi/${subj}_start_left.nii','${subj}_start_left');dtiImportRoiFromNifti('$analysis_folder/ROIs_dwi/${subj}_end_right.nii','${subj}_end_right');dtiImportRoiFromNifti('$analysis_folder/ROIs_dwi/${subj}_end_left.nii','${subj}_end_left')" + +cd $cur_path +} + + +Evaluate_acpc (){ + +cur_path=$(pwd) +subj=$1 + +analysis_folder=/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/$subj/ROIs +cd $analysis_folder + +mrconvert $analysis_folder/$subj\_sr.mif $analysis_folder/$subj\_start_right.nii +mrconvert $analysis_folder/$subj\_sl.mif $analysis_folder/$subj\_start_left.nii +mrconvert $analysis_folder/$subj\_er.mif $analysis_folder/$subj\_end_right.nii +mrconvert $analysis_folder/$subj\_el.mif $analysis_folder/$subj\_end_left.nii + +matlab -r "addpath(genpath('/home/auguser2016/Software/vistasoft-master'));addpath(genpath('/home/auguser2016/Software/mba'));addpath(genpath('/home/auguser2016/Software/encode'));dtiImportRoiFromNifti('$analysis_folder/${subj}_start_right.nii','${subj}_start_right');dtiImportRoiFromNifti('$analysis_folder/${subj}_start_left.nii','${subj}_start_left');dtiImportRoiFromNifti('$analysis_folder/${subj}_end_right.nii','${subj}_end_right');dtiImportRoiFromNifti('$analysis_folder/${subj}_end_left.nii','${subj}_end_left')" + +cd $cur_path +} + + +Evaluate_acpc $1 diff --git a/07_results_evaluation.txt~ b/07_results_evaluation.txt~ new file mode 100644 index 0000000..b6f9509 --- /dev/null +++ b/07_results_evaluation.txt~ @@ -0,0 +1,39 @@ +Evaluate_anatomy (){ + +cur_path=$(pwd) +subj=$1 + +analysis_folder=/home/auguser2016/Projects/Chiasm_RJP_FP/$subj +mkdir $analysis_folder/ROIs +cd $analysis_folder/ROIs + +mrconvert $analysis_folder/$subj\_start_right.mif $analysis_folder/ROIs/$subj\_start_right.nii +mrconvert $analysis_folder/$subj\_start_left.mif $analysis_folder/ROIs/$subj\_start_left.nii +mrconvert $analysis_folder/$subj\_end_right.mif $analysis_folder/ROIs/$subj\_end_right.nii +mrconvert $analysis_folder/$subj\_end_left.mif $analysis_folder/ROIs/$subj\_end_left.nii + +matlab -r "addpath(genpath('/home/auguser2016/Software/vistasoft-master'));addpath(genpath('/home/auguser2016/Software/mba'));addpath(genpath('/home/auguser2016/Software/encode'));dtiImportRoiFromNifti('$analysis_folder/ROIs/${subj}_start_right.nii','${subj}_start_right');dtiImportRoiFromNifti('$analysis_folder/ROIs/${subj}_start_left.nii','${subj}_start_left');dtiImportRoiFromNifti('$analysis_folder/ROIs/${subj}_end_right.nii','${subj}_end_right');dtiImportRoiFromNifti('$analysis_folder/ROIs/${subj}_end_left.nii','${subj}_end_left')" + +cd $cur_path +} + +Evaluate_dwi (){ + +cur_path=$(pwd) +subj=$1 + +analysis_folder=/home/auguser2016/Projects/Chiasm_RJP_FP/$subj +mkdir $analysis_folder/ROIs_dwi +cd $analysis_folder/ROIs_dwi + +mrconvert $analysis_folder/$subj\_sr.mif $analysis_folder/ROIs_dwi/$subj\_start_right.nii +mrconvert $analysis_folder/$subj\_sl.mif $analysis_folder/ROIs_dwi/$subj\_start_left.nii +mrconvert $analysis_folder/$subj\_er.mif $analysis_folder/ROIs_dwi/$subj\_end_right.nii +mrconvert $analysis_folder/$subj\_el.mif $analysis_folder/ROIs_dwi/$subj\_end_left.nii + +matlab -r "addpath(genpath('/home/auguser2016/Software/vistasoft-master'));addpath(genpath('/home/auguser2016/Software/mba'));addpath(genpath('/home/auguser2016/Software/encode'));dtiImportRoiFromNifti('$analysis_folder/ROIs_dwi/${subj}_start_right.nii','${subj}_start_right');dtiImportRoiFromNifti('$analysis_folder/ROIs_dwi/${subj}_start_left.nii','${subj}_start_left');dtiImportRoiFromNifti('$analysis_folder/ROIs_dwi/${subj}_end_right.nii','${subj}_end_right');dtiImportRoiFromNifti('$analysis_folder/ROIs_dwi/${subj}_end_left.nii','${subj}_end_left')" + +cd $cur_path +} + +Evaluate_anatomy $1 diff --git a/08_psychophysics.sh b/08_psychophysics.sh new file mode 100644 index 0000000..0d606fe --- /dev/null +++ b/08_psychophysics.sh @@ -0,0 +1,19 @@ +generate_images { + +subj=$1 + +#Open T1 generate_images + +# Generate random rotation matrix slightly tilting generated images (3 degrees of freedom) + +# Apply matrix to data + +# Load data into mrtrix + +# Select at random one point (z coordinate of intersecting plane) + +# Capture images + +# Repeat + +} diff --git a/Copy_of_06_life_ACPC.m b/Copy_of_06_life_ACPC.m new file mode 100644 index 0000000..9771a27 --- /dev/null +++ b/Copy_of_06_life_ACPC.m @@ -0,0 +1,79 @@ +% Add required folders to path +addpath(genpath('/home/auguser2016/Software/spm12')) +addpath(genpath('/home/auguser2016/Software/vistasoft-master')) +addpath(genpath('/home/auguser2016/Software/mba')) +addpath(genpath('/home/auguser2016/Software/encode')) + +% Define participant +subj ='fe21' + +% Define relevant folders +general_data=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/',subj,'/') +specific_data=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',subj,'/') +current = pwd + +cd(general_data) + +% Build the file names for the diffusion data, the anatomical MRI. +dwiFile = fullfile(strcat(specific_data,subj,'_clean_150mm_aligned_trilin_noMEC.nii.gz')); +dwiFileRepeat = fullfile(strcat(specific_data,subj,'_clean_150mm_aligned_trilin_noMEC.nii.gz')); +t1File = fullfile(strcat(specific_data,subj,'_t1_acpc.nii.gz')); +fgFileName = fullfile(strcat(specific_data,'Tractography/',subj,'_ACT.tck')); +%fgFileName = fullfile(strcat(specific_data,'Tractography/',subj,'_good.tck')); + +% The final connectome and data astructure will be saved with this name: +feFileName = 'test_ensemble_tracking'; + +L = 360; % Discretization parameter +fe = feConnectomeInit(dwiFile,fgFileName,feFileName,[],dwiFileRepeat,t1File,L,[0,1]); + +Niter = 500; +fe = feSet(fe,'fit',feFitModel(feGet(fe,'model'),feGet(fe,'dsigdemeaned'),'bbnnls',Niter,'preconditioner')); + +% Remove non zero weighted fascicles +fg = feGet(fe,'fibers acpc'); +w = feGet(fe,'fiber weights'); +positive_w = w > 0; +fg = fgExtract(fg, positive_w, 'keep'); + +cd(specific_data) +fgWrite(fg,strcat(subj,'_acpc_filtered.mat')) + +% Use MBA to visualize +colors = {[.75 .25 .1]}; +viewCoords = [90,0]; +proportion_to_show = .05; +threshold_length = 10; +slices = {[18 0 0],[0 -40 0],[0 0 -14]}; % Values in mm from AC + +% Prepare the plot of tract +fg = rmfield(fg,'coordspace'); + +% Pick a percentage of fascicles to display (the PN can be too dense for visualization). +%fibs_indx = RandSample(1:length(fg.fibers),round(length(fg.fibers)*proportion_to_show)); +%fg.fibers = fg.fibers(1:1000); + +% Display fascicles and anatomy +fig_h = figure('name','Whole Brain','color','k'); +hold on + +% display anatomy +%t1 = niftiRead(t1File); +%h = mbaDisplayBrainSlice(t1, slices{1}); +%h = mbaDisplayBrainSlice(t1, slices{2}); + +% Display fasciles +set(gca,'visible','off','ylim',[-108 69],'xlim',[-75 75],'zlim',[-45 78],'Color','w') +[~, light_h] = mbaDisplayConnectome(fg.fibers,fig_h,colors{1},'single'); +delete(light_h) +view(viewCoords(1),viewCoords(2)) +light_h = camlight('right'); +lighting phong; +%set(fig_h,'Units','normalized', 'Position',[0.5 .2 .4 .8]); +set(gcf,'Color',[1 1 1]) +drawnow + +%quit + +% save figure to disk +%print(fig_h,'my_figure','-dpng','-r600');%'-dpsc2') \ No newline at end of file diff --git a/old_03_t1_and_dwi_to_acpc.m b/old_03_t1_and_dwi_to_acpc.m new file mode 100644 index 0000000..b1db325 --- /dev/null +++ b/old_03_t1_and_dwi_to_acpc.m @@ -0,0 +1,51 @@ +% Add required folders to path +addpath(genpath('/home/auguser2016/Software/spm12')) +addpath(genpath('/home/auguser2016/Software/vistasoft-master')) +addpath(genpath('/home/auguser2016/Software/mba')) +addpath(genpath('/home/auguser2016/Software/encode')) + +% ACPC registration for T1 in each data set performed manually +% Done manually + +% Perform Diffusion Imaging Processing for all data sets +cd /home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/ +subjects=dir + +for i=5:size(subjects,1) + + id=subjects(i).name + %id='fe21' + % Prepare data in form acceptable by dtiInit + + % ACPC registerred T1 image + anatomy=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/',id,'_t1_acpc.nii.gz') + + % Diffusion Data (motion, eddy and geometrical distortions corrected) + dwi=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/',id,'/',id,'_clean_075mm.nii.gz') + %dwi=strcat('/media/auguser2016/Volume/Test/lw37.nii.gz') + + % Bvecs and bvals obtained from corrected DW images + bvecs=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/',id,'/',id,'_bvecs.bvecs') + bvals=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/',id,'/',id,'_bvals.bvals') + + % Output directory + output=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/') + + % Set up dtiInit paramters specific for our study + dwParams = dtiInitParams; + dwParams.rotateBvecsWithRx = 1; + dwParams.rotateBvecsWithCanXform = 0; + dwParams.bvecsFile = fullfile(bvecs); + dwParams.bvalsFile = fullfile(bvals); + dwParams.eddyCorrect = -1; + dwParams.dwOutMm = [0.75 0.75 0.75]; % change this to 1.5 and process 0.75 anyway + %dwParams.dwOutMm = [1.5 1.5 1.5]; + dwParams.phaseEncodeDir = 2; + dwParams.outDir=output + %dwParams.outDir = fullfile(pwd,'dtiInitPreprocessed'); + + % Start Diffusion Imaging preprocessing. + [dt6FileName, outBaseDir] = dtiInit(dwi,... + anatomy, ... + dwParams); +end diff --git a/option_01_check_preprocessing.sh b/option_01_check_preprocessing.sh new file mode 100644 index 0000000..8b717c0 --- /dev/null +++ b/option_01_check_preprocessing.sh @@ -0,0 +1,81 @@ +Check_preprocessing () { + +subj=$1 + +low_res=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/$subj +high_res=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj + +# View denoising maps. The less anatomical structures are visible the better +mrview $low_res/Temporal/Distortions/$subj\_series* + +# Compare DWI before and after denoising for each encoding direction +mrview $low_res/Temporal/*AP.mif $low_res/Temporal/*AP_denoised.mif +mrview $low_res/Temporal/*PA.mif $low_res/Temporal/*PA_denoised.mif + +# Compare mrtrix gradient for corrected DWI with original gradient table +mrinfo $low_res/Temporal/$subj\_*AP_denoised.mif -export_grad_fsl $low_res/Temporal/bvecs $low_res/Temporal/bvals +gedit $low_res/Temporal/bvecs & +gedit $low_res/Temporal/Distortions/$subj\_gradients & +gedit /home/auguser2016/dMRI_DATA/RAW_DATA/$subj*/s*/*ep2d*/*.bvec +rm $low_res/Temporal/bvecs $low_res/Temporal/bvals + +# Compare DWI before and after dwipreproc +mrview $low_res/Temporal/*AP.mif $low_res/Temporal/*PA.mif $low_res/Temporal/$subj\_denoised_preproc.mif + +# View mask used for bias field correction +mrview $low_res/Temporal/$subj\_denoised_preproc.mif -overlay.load $low_res/$subj\_clean_150mm_mask.mif -overlay.opacity 0.2 -overlay.colourmap 3 + +# View bias field +mrview $low_res/Temporal/Distortions/$subj\_bias_field.mif + +# View DWI before and after bias correction +mrview $low_res/Temporal/$subj\_denoised_preproc.mif $low_res/$subj\_clean_150mm.mif + +# View normal and upsampled images +mrview $low_res/$subj\_clean_150mm.mif $high_res/$subj\_clean_075mm.mif + +# View upsampled mask +mrview $high_res/$subj\_clean_075mm.mif -overlay.load $high_res/$subj\_clean_075mm_mask.mif -overlay.opacity 0.2 -overlay.colourmap 3 + +# View SFR estimation voxels +mrview $high_res/$subj\_clean_075mm.mif -overlay.load $high_res/Temporal/$subj\_clean_075mm_SFR_voxels.mif + +# View SFR +shview $high_res/$subj\_clean_075mm_SFR.txt + +# View anatomical image overlaid on DWI +mrview $high_res/$subj\_clean_075mm_FOD.mif $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz + +# View results of bet2 +mrview $high_res/$subj\_anatomical.nii.gz -overlay.load $high_res/Temporal/$subj\_anatomical_brain.nii.gz -overlay.opacity 0.3 + +# View segmentation +mrview $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.load $high_res/$subj\_075mm_5tt.nii.gz -overlay.opacity 0.3 + +# Correct WM segmentation in chiasm - modify opened ROI and close mrview. Modified WM segmentation image will be shown - repeat previous steps until satisfactory result is obtained. Terminate console to finish process +mrview $high_res/$subj\_075mm_5tt.nii.gz -colourmap 3 -overlay.load $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $high_res/$subj\_5tt_patch.mif + +5ttedit $high_res/$subj\_075mm_5tt.nii.gz -wm $high_res/$subj\_5tt_patch.mif $high_res/$subj\_5tt_for_act_modified.nii.gz -force + +while true; do + mrview $high_res/$subj\_5tt_for_act_modified.nii.gz -colourmap 3 -overlay.load $high_res/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $high_res/$subj\_5tt_patch.mif + 5ttedit $high_res/$subj\_075mm_5tt.nii.gz -wm $high_res/$subj\_5tt_patch.mif $high_res/$subj\_5tt_for_act_modified.nii.gz -force +done + +#data=/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/0.75mm_iso/$subj + +### Correct WM segmentation in chiasm - modify opened ROI and close mrview. Modified WM segmentation image will be shown - repeat previous steps until satisfactory result is obtained. Terminate console to finish process +#mrview $data/$subj\_075mm_5tt.nii.gz -colourmap 3 -overlay.load $data/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $data/$subj\_5tt_patch.mif + +#5ttedit $data/$subj\_075mm_5tt.nii.gz -wm $data/$subj\_5tt_patch.mif $data/$subj\_075mm_5tt_modify.nii.gz -force + +#while true; do +# mrview $data/$subj\_075mm_5tt_modify.nii.gz -colourmap 3 -overlay.load $data/$subj\_075mm_anatomical_coreg2nodif_undist_2.nii.gz -overlay.colourmap 0 -overlay.opacity 0.8 -roi.load $data/$subj\_5tt_patch.mif +# 5ttedit $data/$subj\_075mm_5tt.nii.gz -wm $data/$subj\_5tt_patch.mif $data/$subj\_075mm_5tt_modify.nii.gz -force +#done + +} + +#for i in /home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/*; do Check_preprocessing ${i: -4}; done + +Check_preprocessing $1 diff --git a/option_03_working_on_original_lw.m b/option_03_working_on_original_lw.m new file mode 100644 index 0000000..6ae4082 --- /dev/null +++ b/option_03_working_on_original_lw.m @@ -0,0 +1,49 @@ +% Add required folders to path +addpath(genpath('/home/auguser2016/Software/spm12')) +addpath(genpath('/home/auguser2016/Software/vistasoft-master')) +addpath(genpath('/home/auguser2016/Software/mba')) +addpath(genpath('/home/auguser2016/Software/encode')) + +% ACPC registration for T1 in each data set performed manually +% Done manually + +% Perform Diffusion Imaging Processing for all data sets +cd /home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/ +subjects=dir + +for i=5:size(subjects,1) + + id=subjects(i).name + % Prepare data in form acceptable by dtiInit + + % ACPC registerred T1 image + anatomy=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/',id,'_t1_acpc.nii.gz') + + % Diffusion Data (motion, eddy and geometrical distortions corrected) + dwi=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_clean_150mm.nii.gz') + + % Bvecs and bvals obtained from corrected DW images + bvecs=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_bvecs.bvecs') + bvals=strcat('/home/auguser2016/dMRI_DATA/PREPROCESSED_DATA/1.5mm_iso/',id,'/',id,'_bvals.bvals') + + % Output directory + output=strcat('/home/auguser2016/Projects/Chiasm_RJP_FP_ACPC/',id,'/') + + % Set up dtiInit paramters specific for our study + dwParams = dtiInitParams; + dwParams.rotateBvecsWithRx = 1; + dwParams.rotateBvecsWithCanXform = 0; + dwParams.bvecsFile = fullfile(bvecs); + dwParams.bvalsFile = fullfile(bvals); + dwParams.eddyCorrect = -1; + dwParams.dwOutMm = [0.75 0.75 0.75]; % change this to 1.5 and process 0.75 anyway + %dwParams.dwOutMm = [1.5 1.5 1.5]; + dwParams.phaseEncodeDir = 2; + dwParams.outDir=output + %dwParams.outDir = fullfile(pwd,'dtiInitPreprocessed'); + + % Start Diffusion Imaging preprocessing. + [dt6FileName, outBaseDir] = dtiInit(dwi,... + anatomy, ... + dwParams); +end