dHCP_neo_dMRI_importFiles.sh 7.24 KB
Newer Older
Matteo Bastiani's avatar
Matteo Bastiani committed
1 2 3 4 5 6 7 8 9
#!/bin/bash


set -e
echo -e "\n START: Importing files"


if [ "$6" == "" ];then
    echo ""
10
    echo "usage: $0 <Data folder> <Data filename> <Output folder> <Acquisition protocol> <Number of volumes> <Number of b0s>"
Matteo Bastiani's avatar
Matteo Bastiani committed
11 12
    echo "       Data folder: Path to the data file"
    echo "       Data filename: File name for the raw data"
13
    echo "       Output folder: Path to the subject output folder"
Matteo Bastiani's avatar
Matteo Bastiani committed
14 15 16 17 18 19 20 21 22 23
    echo "       Acqusition protocol: Text file with acquisition parameters"
    echo "       Number of volumes: Number of acquired volumes"
    echo "       Number of b0s: Number of b0s from each phase encoding block to use when estimating the fieldmap"
    echo ""
    exit 1
fi


dataFolder=$1      # Path to data file
dataFile=$2        # Data file name
24
outFolder=$3       # Subject output folder
Matteo Bastiani's avatar
Matteo Bastiani committed
25 26 27 28
acqProt=$4         # Acquisition protcol
nVols=$5           # Number of acquired volumes
noB0s=$6           # Number of B0 volumes for each PE direction used to estimate distortions with TOPUP

29 30 31 32
prepFolder=${outFolder}/PreProcessed
rawFolder=${outFolder}/raw
mkdir -p ${rawFolder}/tmpData

Matteo Bastiani's avatar
Matteo Bastiani committed
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

#============================================================================
# Separate b0 volumes according to their PE direction and obtain 
# bvals and bvecs files. Round bvals to identify shells.
#============================================================================
idxLR=0
idxRL=0
idxAP=0
idxPA=0
LRvols=-1
RLvols=-1
APvols=-1
PAvols=-1
no_LR=0
no_RL=0
no_AP=0
no_PA=0
idxb400=0
idxb1000=0
idxb2600=0
53 54
bshells=(100000)
n_b=0
Matteo Bastiani's avatar
Matteo Bastiani committed
55 56
i=0

57
echo -n > ${rawFolder}/tmpData/eddyIndex.txt
Matteo Bastiani's avatar
Matteo Bastiani committed
58 59 60 61 62 63 64
while read line; do
    bvec_x[i]=`echo $line | awk {'print $1'}`
    bvec_y[i]=`echo $line | awk {'print $2'}`
    bvec_z[i]=`echo $line | awk {'print $3'}`
    bval[i]=`echo $line | awk {'print $4'}`
    pedir[i]=`echo $line | awk {'print $5'}`
    rotime[i]=`echo $line | awk {'print $6'}`
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
    echo -n "${pedir[i]} " >> ${rawFolder}/tmpData/eddyIndex.txt
    if [ ${bval[i]} -lt 100 ]; then      #b0
        first_b0=${i}
        if [ ${pedir[i]} -eq 1 ]; then   #LR
            LRvols[idxLR]=$i
            idxLR=$(($idxLR + 1))	
        elif [ ${pedir[i]} -eq 2 ]; then #RL
            RLvols[idxRL]=$i
            idxRL=$(($idxRL + 1))
        elif [ ${pedir[i]} -eq 3 ]; then #AP
            APvols[idxAP]=$i
            idxAP=$(($idxAP + 1))
        elif [ ${pedir[i]} -eq 4 ]; then #PA
            PAvols[idxPA]=$i
            idxPA=$(($idxPA + 1))
        fi
    else                                 #dw
        if [ ${pedir[i]} -eq 1 ]; then   #LR
            no_LR=$(($no_LR + 1))
        elif [ ${pedir[i]} -eq 2 ]; then #RL
            no_RL=$(($no_RL + 1))
        elif [ ${pedir[i]} -eq 3 ]; then #AP
            no_AP=$(($no_AP + 1))
        elif [ ${pedir[i]} -eq 4 ]; then #PA
            no_PA=$(($no_PA + 1))
        fi
    fi
    # Identify unique shells
    b_flag=0
    for ub in "${bshells[@]}"; do 
        j=`echo ${bval[i]} | awk -F"E" 'BEGIN{OFMT="%10.10f"} {print $1 * (10 ^ $2)}' `  # Round bval
	    j=${j/.*}   
	    b_diff=`echo "${j} - ${ub}" | bc | awk ' { if($1>=0) { print $1} else {print $1*-1 }}'`
        if [ ${b_diff} -lt 100 ]; then
            b_flag=1
            break
        fi
    done
    if [ "${b_flag}" == 0 ]; then
        bshells[n_b]=${bval[i]}
        n_b=$((${n_b} + 1))
Matteo Bastiani's avatar
Matteo Bastiani committed
106
    fi
107

Matteo Bastiani's avatar
Matteo Bastiani committed
108 109
    i=$(($i + 1))
    if [ $i -eq ${nVols} ]; then
110
	    break
Matteo Bastiani's avatar
Matteo Bastiani committed
111 112 113
    fi
done < ${acqProt}

114 115 116 117 118 119
bshells=($(echo "${bshells[@]}" | tr ' ' '\n' | sort -n -u | tr '\n' ' '))
echo "${bshells[@]}" > ${rawFolder}/shells
echo "Found ${n_b} shells: ${bshells[@]} s/mm^2"
un_pedirs=(`echo "${pedir[@]}" | tr ' ' '\n' | sort -n -u | tr '\n' ' '`)
echo "${un_pedirs[@]}" > ${rawFolder}/pedirs

Matteo Bastiani's avatar
Matteo Bastiani committed
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136

#============================================================================
# Write json file for QC
#============================================================================
echo "{" > ${prepFolder}/dataImport.json

echo "   \"no_LR_vols\": $no_LR," >> ${prepFolder}/dataImport.json
echo "   \"no_RL_vols\": $no_RL," >> ${prepFolder}/dataImport.json
echo "   \"no_AP_vols\": $no_AP," >> ${prepFolder}/dataImport.json
echo "   \"no_PA_vols\": $no_PA" >> ${prepFolder}/dataImport.json

echo "}" >> ${prepFolder}/dataImport.json


#============================================================================
# Write bvecs, bvals and eddy acquisition parameters file.
#============================================================================
137 138 139 140
echo "${bvec_x[@]}" > ${rawFolder}/tmpData/orig_bvecs
echo "${bvec_y[@]}" >> ${rawFolder}/tmpData/orig_bvecs
echo "${bvec_z[@]}" >> ${rawFolder}/tmpData/orig_bvecs
echo "${bval[@]}" > ${rawFolder}/tmpData/bvals
Matteo Bastiani's avatar
Matteo Bastiani committed
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162

if [ $LRvols -ne -1 ]; then
    echo -1 0 0 ${rotime[${LRvols[0]}]} > ${prepFolder}/eddy/acqparamsUnwarp.txt
else
    echo -1 0 0 0.05 > ${prepFolder}/eddy/acqparamsUnwarp.txt
fi
if [ $RLvols -ne -1 ]; then
    echo 1 0 0 ${rotime[${RLvols[0]}]} >> ${prepFolder}/eddy/acqparamsUnwarp.txt
else
    echo 1 0 0 0.05 >> ${prepFolder}/eddy/acqparamsUnwarp.txt
fi
if [ $APvols -ne -1 ]; then
    echo 0 -1 0 ${rotime[${APvols[0]}]} >> ${prepFolder}/eddy/acqparamsUnwarp.txt
else
    echo 0 -1 0 0.05 >> ${prepFolder}/eddy/acqparamsUnwarp.txt
fi
if [ $PAvols -ne -1 ]; then
    echo 0 1 0 ${rotime[${PAvols[0]}]} >> ${prepFolder}/eddy/acqparamsUnwarp.txt
else
    echo 0 1 0 0.05 >> ${prepFolder}/eddy/acqparamsUnwarp.txt
fi

163

Matteo Bastiani's avatar
Matteo Bastiani committed
164
#============================================================================
165 166
# Reorient raw data and bvecs to to match the approximate orientation 
# of the standard template images (MNI152)
Matteo Bastiani's avatar
Matteo Bastiani committed
167
#============================================================================
168 169
${FSLDIR}/bin/fslreorient2std ${dataFolder}/${dataFile} ${rawFolder}/tmpData/data
${FSLDIR}/bin/fslreorient2std ${dataFolder}/${dataFile} > ${rawFolder}/tmpData/raw2std.mat
Matteo Bastiani's avatar
Matteo Bastiani committed
170

171
${scriptsFolder}/utils/rotateBvecs.sh ${rawFolder}/tmpData/orig_bvecs ${rawFolder}/tmpData/raw2std.mat ${rawFolder}/tmpData/bvecs
Matteo Bastiani's avatar
Matteo Bastiani committed
172 173 174


#============================================================================
175 176
# If more than 1 PE direction has been acquired, Identify best b0s for each 
# PE direction; otherwise, set the first b0 as the reference volume.
Matteo Bastiani's avatar
Matteo Bastiani committed
177
#============================================================================
178 179 180 181 182 183 184
unique_pedirs=(`cat ${rawFolder}/pedirs`)
if [ `echo ${#unique_pedirs[@]}` -gt 1 ]; then
    echo "More than 1 phase encoding direction detected. Selecting best b0 volumes for topup"
    ${scriptsFolder}/utils/pickBestB0s ${rawFolder}/tmpData/data ${rawFolder}/tmpData/bvals ${rawFolder}/tmpData/eddyIndex.txt ${rotime[0]} ${noB0s} ${prepFolder}/topup
else
    echo "1 phase encoding direction detected. Setting first b0 as reference volume"
    echo "${first_b0}" > ${prepFolder}/topup/ref_b0.txt
Matteo Bastiani's avatar
Matteo Bastiani committed
185 186 187 188
fi


#============================================================================
189
# Sort raw data based on the selected reference volume
Matteo Bastiani's avatar
Matteo Bastiani committed
190
#============================================================================
191
${scriptsFolder}/utils/sortData ${rawFolder}/tmpData/data `cat ${prepFolder}/topup/ref_b0.txt` ${rawFolder} ${rawFolder}/tmpData/bvals ${rawFolder}/tmpData/bvecs ${rawFolder}/tmpData/eddyIndex.txt
Matteo Bastiani's avatar
Matteo Bastiani committed
192 193 194


#============================================================================
195
# Clean unnecessary files
Matteo Bastiani's avatar
Matteo Bastiani committed
196
#============================================================================
197
rm -rf ${rawFolder}/tmpData
Matteo Bastiani's avatar
Matteo Bastiani committed
198 199 200 201


echo -e "\n END: Importing files"