Skip to content
Snippets Groups Projects
Commit 0bfa8784 authored by Rita Nunes's avatar Rita Nunes
Browse files

Initial import

parents
No related branches found
Tags misc_scripts-1_0
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Makefile 0 → 100644
+ 7
− 0
View file @ 0bfa8784
include ${FSLCONFDIR}/default.mk
PROJNAME = misc_scripts
SCRIPTS = regscript remove_vols replace_and_average_fmrib
all:
regscript 0 → 100755
+ 352
− 0
View file @ 0bfa8784
#!/usr/local/bin/bash
# preprocessing for DTI images acquired at FMRIB
####################################################################
#This script performs registration of all the volumes
#in a given dti dataset to a reference non-diffusion-weighted image
#(scaling, translation along the PE direction and shearing)
#Brain Extraction is also performed
#with a fractional intensity treshold of 0.025
####################################################################
Usage() {
echo "REGSCRIPT - Preprocessing for DTI images acquired at FMRIB "
echo ""
echo "USAGE: regscript <username> <study_number> <output_dir> [number_acqs] [type_of_correction] [series_numb_dti] [zero_pad] [study_dir]"
echo ""
echo "default number_acqs is 1"
echo ""
echo "The FID files should normally be in"
echo " ~username/scratch/MRdata/study_number/"
echo ""
echo "type_of_correction (optional):"
echo " mc - correct for eddy currents and motion (default)"
echo " bvecs will only be rotated if number_acqs is 1"
echo " as the amount of rotation corresponding to a given"
echo " volume for different acquisitions would not be the same"
echo " and it wouldn't be possible to simply average the data"
echo " ed - correct only for eddy currents (3 degrees of freedom)"
echo ""
echo "series_numb_dti (optional)- default is 2"
echo " e.g. series_2.1 - non-DW image (nodif)"
echo " series_2.2 - DTI (acquisition 1)"
echo " series_2.3 - DTI (acquisition 2) ..."
echo ""
echo "zero_pad (optional) - default is off"
echo " set it on to run load_varian"
echo " with the option -hkspad "
echo " you might wish to try this if the"
echo " half k-space reconstruction has failed"
echo ""
echo "Indicate study_dir if the files are in a different directory"
echo ""
echo "OUTPUTS"
echo "-----------------------------------------------------------------------"
echo " non-DWI: nodif"
echo " DTI after registration: data (if number_acqs equal to 1)"
echo " data1, data2... (several repetitions)"
echo " Brain mask: nodif_brain_mask"
echo " T1-structural (optional - only if series_numb_dti is 2): struct"
echo " T2-structural (optional - only if series_numb_dti is 2): T2-weighted"
echo " Diffusion directions: bvals and bvecs"
echo ""
exit
}
[ "$1" = "" ] && Usage
[ "$2" = "" ] && Usage
[ "$3" = "" ] && Usage
username=$1
study_number=$2
output_dir=`pwd`/$3
if [ $# -ge 4 ]; then
number_acqs=$4
echo Assuming number_acqs is $number_acqs
else
number_acqs=1;
echo Assuming number_acqs is 1
fi
if [ $# -ge 5 ]; then
type_cor=$5
else
type_cor='mc'
fi
if [ $# -ge 6 ]; then
series_numb=$6
echo series_numb_dti is $series_numb
else
series_numb=2;
echo series_numb_dti is $series_numb
fi
option='-hks'
if [ $# -ge 7 ]; then
check=$7
check_option=`echo $check`
case $check_option
in
on) option='-hkspad';;
off) option='-hks';;
*) option='-hks';;
esac
fi
echo load_varian will run with option $option
if [ $# -ge 8 ]; then
study_dir=$8
else
study_dir=/usr/people/${username}/scratch/MRdata
fi
nodiffid=${study_dir}/${study_number}/series_${series_numb}.1
dtifid=${study_dir}/${study_number}/series_${series_numb}.
structfid=${study_dir}/${study_number}/series_3.1
T2fid=${study_dir}/${study_number}/series_4.1
if [ ! -e $output_dir ]; then
mkdir $output_dir
fi
cd $output_dir
echo output_dir is $output_dir
cp ~dtiuser/etc/dof3 .
#Load FID files
if [ -e ${nodiffid}.fid ]; then
echo Loading nodif
load_varian $nodiffid $output_dir/nodif -auto -float -buo -egrv $option
else
echo $nodiffid.fid: FID not found
exit;
fi
nvols_nodif=`avwval $output_dir/nodif dim4`;
if [ $nvols_nodif -gt 1 ]; then
echo nodif has $nvols_nodif volumes - file name will be changed to data
echo first volume will be saved as nodif
mv nodif.hdr dti1.hdr
mv nodif.img dti1.img
avwsplit dti1
mv vol0000.hdr nodif.hdr
mv vol0000.img nodif.img
rm vol* -f
count=1
image_numb=2
scan_numb=2
else
count=0
image_numb=1
scan_numb=2
fi
while [ $count -lt $number_acqs ]; do
if [ -e ${dtifid}${scan_numb}.fid ]; then
echo Loading dti${image_numb}
load_varian ${dtifid}${scan_numb}.fid ${output_dir}/dti${image_numb} -auto -float -buo -egrv $option
else
echo ${dtifid}${image_numb}.fid: file not found
exit;
fi
count=`expr $count + 1`;
image_numb=`expr $image_numb + 1`;
scan_numb=`expr $scan_numb + 1`;
done
scan_numb=`expr $scan_numb - 1`;
if [ $series_numb -eq 2 ]; then
if [ -e ${structfid}.fid ]; then
echo Loading structural
load_varian ${structfid} ${output_dir}/struct -auto -float
fi
if [ -e ${T2fid}.fid ]; then
echo Loading T2-weighted
load_varian $T2fid ${output_dir}/T2_weighted -auto -float
fi
fi
#Registers each volume in dti to the reference image nodif
nvols=`avwval dti1 dim4`;
echo $nvols gradient directions
count_acqs=0
correc=`echo $type_cor`;
case $correc
in
mc) echo Correcting for eddy currents and motion
if [ ! -e matrices ]; then
mkdir matrices
fi
while [ $count_acqs -lt $number_acqs ]; do
tmp_acqs=`expr $count_acqs + 1`;
avwsplit dti${tmp_acqs}
if [ $nvols -lt 9 ]; then
count=0
while [ ${count} != ${nvols} ];do
echo FLIRT on volume ${count}
flirt -in vol000$count -ref nodif -out rvol000$count -nosearch -omat matrices/mat
avscale matrices/mat vol000$count >matrices/tmp
head -4 matrices/tmp | tail -3 > matrices/MAT000${count}${tmp_acqs}
count=`expr $count + 1`;
done
else
count=0
while [ $count != 10 ];do
echo FLIRT on volume ${count}
flirt -in vol000$count -ref nodif -out rvol000$count -nosearch -omat matrices/mat
avscale matrices/mat vol000$count >matrices/tmp
head -4 matrices/tmp | tail -3 > matrices/MAT000${count}${tmp_acqs}
count=`expr $count + 1`;
done
while [ ${count} != ${nvols} ];do
echo FLIRT on volume ${count}
flirt -in vol00${count} -ref nodif -out rvol00${count} -nosearch -omat matrices/mat
avscale matrices/mat vol00$count >matrices/tmp
head -4 matrices/tmp | tail -3 > matrices/MAT00${count}${tmp_acqs}
rm -f matrices/mat
rm -f matrices/tmp
count=`expr $count + 1`;
done
fi
#Clean up unecessary files
rm vol* -f
rm dti${tmp_acqs}* -f
if [ $number_acqs -eq 1 ]
then
echo creating file: data
avwmerge -t data rvol00*.hdr
else
echo creating file: data${tmp_acqs}
avwmerge -t data${tmp_acqs} rvol00*.hdr
fi
#Clean up unecessary files
rm rvol* -f
count_acqs=`expr $count_acqs + 1`;
done
rm -f matrices/mat;
rm -f matrices/tmp;;
ed) echo Correcting only for eddy currents
while [ $count_acqs -lt $number_acqs ]; do
tmp_acqs=`expr $count_acqs + 1`;
avwsplit dti${tmp_acqs}
if [ $nvols -lt 9 ]; then
count=0
while [ ${count} != ${nvols} ];do
echo FLIRT on volume ${count}
flirt -schedule dof3 -in vol000$count -ref nodif -out rvol000$count -nosearch
count=`expr $count + 1`;
done
else
count=0
while [ $count != 10 ];do
echo FLIRT on volume ${count}
flirt -schedule dof3 -in vol000$count -ref nodif -out rvol000$count -nosearch
count=`expr $count + 1`;
done
while [ ${count} != ${nvols} ];do
echo FLIRT on volume ${count}
flirt -schedule dof3 -in vol00${count} -ref nodif -out rvol00${count} -nosearch
count=`expr $count + 1`;
done
fi
#Clean up unecessary files
rm vol* -f
rm dti${tmp_acqs}* -f
if [ $number_acqs -eq 1 ]
then
echo creating file: data
avwmerge -t data rvol00*.hdr
else
echo creating file: data${tmp_acqs}
avwmerge -t data${tmp_acqs} rvol00*.hdr
fi
#Clean up unecessary files
rm rvol* -f
count_acqs=`expr $count_acqs + 1`;
done;;
*) echo Type_of_correction must be either mc or ed
exit;;
esac
#Brain Extraction
echo BE
bet nodif mask -f 0.025 -m
rm mask.hdr -f
rm mask.img -f
mv mask_mask.hdr nodif_brain_mask.hdr
mv mask_mask.img nodif_brain_mask.img
rm dof3 -f
#Copy relevant bvals and bvecs file
n_vols=`echo $nvols`;
case $n_vols
in
7) cp ~dtiuser/etc/FMRIB_bvals_bvecs/bvals_6 bvals;
cp ~dtiuser/etc/FMRIB_bvals_bvecs/bvecs_6 bvecs;
filename='gdata7';;
13) cp ~dtiuser/etc/FMRIB_bvals_bvecs/bvals_12 bvals;
cp ~dtiuser/etc/FMRIB_bvals_bvecs/bvecs_12 bvecs;
filename='gdata13';;
63) cp ~dtiuser/etc/FMRIB_bvals_bvecs/bvals_60 bvals;
cp ~dtiuser/etc/FMRIB_bvals_bvecs/bvecs_60 bvecs;
filename='gdata63';;
*) echo 'bvals and bvecs not available for that # of directions';
exit;;
esac
#The bvecs will only be rotated for 1 acquisition
#The amount of rotation corresponding to a given volume for
#different acquisitions would be different and it wouldn't be possible
#to simply average the data
if [ $number_acqs -eq 1 ]; then
echo Opening Matlab
case $correc
in
mc)
echo "addpath('/usr/people/dtiuser/etc/FMRIB_bvals_bvecs/')">tmp1
if [ $nvols_nodif -gt 1 ]; then
echo -n "rotate_bvecs($number_acqs,$n_vols,'$nodiffid')" > tmp2
else
echo -n "rotate_bvecs($number_acqs,$n_vols,'${dtifid}${scan_numb}.fid')" > tmp2;
fi
echo >tmp3
cat tmp1 tmp2 tmp3 >tmp.m
cat tmp.m |matlab;
rm tmp1 tmp2 tmp3 tmp.m -f;;
*) exit;;
esac
fi
echo Done!
remove_vols 0 → 100755
+ 89
− 0
View file @ 0bfa8784
#!/usr/local/bin/bash
# preprocessing for DTI images acquired at FMRIB
Usage() {
echo "REMOVE_VOLS - Remove dodgy volumes from DTI data acquired at FMRIB "
echo ""
echo "USAGE: remove_vols <input_file> <output_file> <dodgyvols> "
echo ""
echo "example: remove_vols data data_rv 4 19 35"
echo ""
echo "Note that the first volume is 0 and not 1"
echo ""
exit
}
[ "$1" = "" ] && Usage
[ "$2" = "" ] && Usage
[ "$3" = "" ] && Usage
input_file=$1
output_file=$2
current_dir=`pwd`
if [ ! -e $input_file.hdr ]; then
echo Cant find $input_file
exit
fi
if [ ! -e bvals ]; then
echo Cant find bvals
exit
fi
if [ ! -e bvecs ]; then
echo Cant find bvecs
exit
fi
avwsplit $input_file.hdr
vol=0
num_vols=`expr $# - 2`
t=""
rm dodgy_vols -f
echo -n $t >dodgy_vols
while [ $vol -lt $num_vols ]; do
echo -n $t>tmp1
echo -n $t>tmp2
v=$3
shift
if [ -e vol000${v}.hdr ]; then
rm vol000${v}.hdr -f
rm vol000${v}.img -f
fi
if [ -e vol00${v}.hdr ]; then
rm vol00${v}.hdr -f
rm vol00${v}.img -f
fi
echo $v > tmp1
cat dodgy_vols tmp1 > tmp2
cat tmp2 > dodgy_vols
vol=`expr $vol + 1`
done
echo "addpath('/usr/people/dtiuser/etc/FMRIB_bvals_bvecs/')">tmp3
echo "remove_vols('$current_dir')">tmp4
cat tmp3 tmp4 >tmp.m
cat tmp.m | matlab
rm tmp* -f
echo "Merging Volumes"
avwmerge -t $output_file vol*.hdr
rm vol* -f
rm dodgy_vols -f
replace_and_average_fmrib 0 → 100755
+ 327
− 0
View file @ 0bfa8784
#!/usr/local/bin/bash
#averages multiple repeats and simultaneously remove dodgy volumes
#to be used after regscript - for DTI images acquired at FMRIB
Usage() {
echo ""
echo "replace_and_average_fmrib - Preprocessing for DTI images acquired at FMRIB "
echo " run after regscript to average multiple repeats"
echo ""
echo "if necessary, indicate dodgy volumes to be removed"
echo ""
echo "_________________________________________________________________________"
echo "Usage: replace_and_average_fmrib <numb_acqs> volno volno ... volno "
echo ""
echo "output file will be called data"
echo "Note that the first volume is 0 and not 1"
echo ""
echo "Example: if you want to remove vols 4 and 5 from data1"
echo " vols 10 and 19 from data2"
echo " vol 40 from data3"
echo " having acquired 3 repeats, 60 directions"
echo ""
echo " to the volumes corresponding to data2 (10 and 19)"
echo " you should add the total number of volumes in each file,"
echo " in this case 63 (= 3 non-DW + 60 direc)"
echo " 126 (= 2*63) to the ones corresponding to data3 and so on..."
echo ""
echo "you should therefore type something like - "
echo " replace_and_average_fmrib 3 4 5 73 82 166"
echo "_________________________________________________________________________"
exit
}
#===================================================================
rem_vols() {
#===================================================================
# Arg_1 = input_file
# Arg_2 = num_acqs
# Arg_3 ... Arg_n dodgy_vols
input_file=$1
numb_acqs=$2
shift 2
dodgy_vols=`cat dodgy_vols`
echo "removing vols $dodgy_vols"
current_dir=`pwd`
acq=1
while [ $acq -le $numb_acqs ];do
if [ ! -e ${input_file}${acq}.hdr ]; then
echo Cant find $input_file
exit
fi
acq=`expr $acq + 1`;
done
if [ ! -e bvals ]; then
echo Cant find bvals
exit
fi
if [ ! -e bvecs ]; then
echo Cant find bvecs
exit
fi
nvols=`avwval ${input_file}1 dim4`;
acq=1
while [ $acq -le $numb_acqs ]; do
echo splitting $input_file$acq
avwsplit $input_file$acq.hdr
volume=0;
if [ $nvols -lt 10 ]; then
while [ $volume -lt $nvols ]; do
mv vol000$volume.hdr vol000${volume}_${acq}.hdr
mv vol000$volume.img vol000${volume}_${acq}.img
volume=`expr $volume + 1`
done
else
volume=0
while [ $volume -lt 10 ]; do
mv vol000${volume}.hdr vol000${volume}_${acq}.hdr
mv vol000${volume}.img vol000${volume}_${acq}.img
volume=`expr $volume + 1`
done
while [ $volume -lt $nvols ]; do
mv vol00${volume}.hdr vol000${volume}_${acq}.hdr
mv vol00${volume}.img vol000${volume}_${acq}.img
volume=`expr $volume + 1`
done
fi
acq=`expr $acq + 1`;
done
vol=0
num_dodgy_vols=`expr $#`
while [ $vol -lt $num_dodgy_vols ]; do
v=$1
shift
if [ -e vol000${v}_1.hdr ]; then
rm vol000${v}_*.hdr -f
rm vol000${v}_*.img -f
fi
if [ -e vol00${v}_1.hdr ]; then
rm vol00${v}_*.hdr -f
rm vol00${v}_*.img -f
fi
vol=`expr $vol + 1`
done
acq=1
while [ $acq -le $numb_acqs ]; do
rm ${input_file}${acq}*
avwmerge -t data$acq vol*_$acq.hdr
acq=`expr $acq + 1`;
done
matlab -nodisplay -nojvm -nosplash 1> matlab.out1 2>&1 <<EOF
addpath('/usr/people/dtiuser/etc/FMRIB_bvals_bvecs/')
remove_vols('$current_dir')
EOF
rm matlab.out1 -f
rm vol0* -f
}
#===================================================================
check_vols() {
#===================================================================
# Arg_1 = n_vols
# Arg_2 = num_acqs
n_vols=$1
numb_acqs=$2
matlab -nodisplay -nojvm -nosplash 1> matlab.out1 2>&1 <<EOF
addpath('/usr/people/dtiuser/etc/FMRIB_bvals_bvecs/')
check_vols($n_vols,$numb_acqs)
EOF
rm matlab.out1 -f
}
#===================================================================
[ "$1" = "" ] && Usage
numb_acqs=$1
if [ $# -gt 1 ]; then
shift 1
if [ $numb_acqs -lt 2 ]; then
echo "replacevols_fmrib should be used when several repeats were acquired"
echo "use remove_vols instead if you want to remove dodgy volumes from single datasets"
fi
echo -n ""> vol_list
for i in $@ ; do
vol=$i
echo -n "$vol " >> vol_list
shift
done
nvols=`avwval data1 dim4`;
n_vols=`echo $nvols`
dodgy_volumes=`cat vol_list`
echo dodgy volumes: $dodgy_volumes
check_vols $n_vols $numb_acqs
if [ -e dodgy_vols ]; then
dodgy=`cat dodgy_vols`
rem_vols data $numb_acqs $dodgy
else
echo No volumes need to be removed
fi
nvols=`avwval data1 dim4`;
echo output will have $nvols volumes
dodgy_volumes=`cat vol_list`
dodgy_volumes=`echo $dodgy_volumes`
if [ -e dodgy_vols ]; then
if [ "$dodgy_volumes" != "" ]; then
echo dogdy volumes updated taking into account removed volumes
echo dodgy volumes are now: $dodgy_volumes
fi
fi
if [ "$dodgy_volumes" != "" ]; then
acq=1
echo "">file_list
while [ $acq -le $numb_acqs ];do
echo data$acq.hdr >> file_list
acq=`expr $acq + 1`;
done
file_list=`cat file_list`
avwmerge -t bigdata $file_list
rm file_list -f
echo "0 $nvols" >average_file
acq=1
first=0
second=$nvols
while [ $acq -lt $numb_acqs ];do
first=`expr $first + $nvols`;
echo "$first $second">>average_file;
acq=`expr $acq + 1`;
done
echo "">>average_file
acq=0
n_vols=`echo $nvols`
if [ $n_vols -lt 63 ]; then
b0=0;
while [ $acq -lt $numb_acqs ];do
echo -n "$b0 " >> average_file
acq=`expr $acq + 1`;
b0=`expr $b0 + $n_vols`;
done;
else
b0=0;b0_1=1;b0_2=2;echo ""
while [ $acq -lt $numb_acqs ]; do
echo -n "$b0 $b0_1 $b0_2 " >> average_file
acq=`expr $acq + 1`;
b0=`expr $b0 + $n_vols`;
b0_1=`expr $b0_1 + $n_vols`;
b0_2=`expr $b0_2 + $n_vols`;
done;
fi
replacevols bigdata average_file $dodgy_volumes tmp_file
if [ -e tmp_file.img.gz ]; then
gunzip tmp_file.img.gz
if [ -e tmp_file.hdr.gz ]; then
gunzip tmp_file.hdr.gz
fi
fi
nvols=`avwval tmp_file dim4`;
n_vols=`echo $nvols`
volsacqs=`avwval data1 dim4`;
vols_per_acqs=`echo $volsacqs`;
rm bigdata* -f
avwroi tmp_file data 0 ${vols_per_acqs};
avwmaths data -mul 0 data;
acq=0;
while [ ${acq} -lt ${numb_acqs} ]; do
vol=0
while [ ${vol} -lt ${vols_per_acqs} ]; do
v_tmp=`expr $vol + ${vols_per_acqs} \* ${acq}`
if [ $vol -lt 10 ]; then
avwroi tmp_file vol${acq}_000${vol} $v_tmp 1;
else
avwroi tmp_file vol${acq}_00${vol} $v_tmp 1;
fi
vol=`expr $vol + 1`;
done
avwmerge -t data${acq}_all vol${acq}_*.hdr
avwmaths data${acq}_all -div ${numb_acqs} data${acq}_all;
add4d data data${acq}_all data;
rm data${acq}_all* vol${acq}* -f
acq=`expr $acq + 1`;
done
rm tmp*
rm *_all* -f
rm average_file -f
else
echo averaging the data
volsacqs=`avwval data1 dim4`;
vols_per_acqs=`echo $volsacqs`;
avwroi data1 data 0 ${vols_per_acqs};
avwmaths data -mul 0 data;
acq=1;
while [ ${acq} -le ${numb_acqs} ]; do
avwmaths data${acq} -div ${numb_acqs} data${acq};
add4d data data${acq} data;
acq=`expr $acq + 1`;
done
fi
rm dodgy_vols -f
rm vol_list -f
rm vol*_* -f
else
echo no dodgy volumes
echo averaging the data
volsacqs=`avwval data1 dim4`;
vols_per_acqs=`echo $volsacqs`;
avwroi data1 data 0 ${vols_per_acqs};
avwmaths data -mul 0 data;
acq=1;
while [ ${acq} -le ${numb_acqs} ]; do
avwmaths data${acq} -div ${numb_acqs} data${acq};
add4d data data${acq} data;
acq=`expr $acq + 1`;
done
fi
echo extracting nodif image from averaged data
avwsplit data
mv vol0000.hdr nodif.hdr -f
mv vol0000.img nodif.img -f
rm tmp* -f
rm vol* -f
rita_uncertainty_gating 0 → 100755
+ 19
− 0
View file @ 0bfa8784
#!/usr/local/bin/bash
main_dir='/home/fs10/dtiuser/RITA/gating'
for direc in 03Aug04_SL
do
cd $direc
for a in 13
do
#avwroi data1 gated_$a 0 128 0 128 $a 1
# avwroi data2 notgated_$a 0 128 0 128 $a 1
diffmodmhgibbs -k notgated_${a} -o gated_$a -m mask_sl$a -r bvecs -b bvals -B 1000 -S 20200 --saverat=20 --scj --noprec;
# diffmodmhgibbs -k notgated_${a} -o notgated_$a -m mask_sl$a -r bvecs -b bvals -B 1000 -S 20200 --saverat=20 --scj;
done
cd $main_dir
done
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment