Codebase list bart / 6e489f9
add script for GRASP Martin Uecker 8 years ago
1 changed file(s) with 230 addition(s) and 0 deletion(s). Raw diff Collapse all Expand all
+230
-0
scripts/grasp.sh less more
0 #!/bin/bash
1 # Copyright 2015. The Regents of the University of California.
2 # All rights reserved. Use of this source code is governed by
3 # a BSD-style license which can be found in the LICENSE file.
4 #
5 # Authors:
6 # 2015 Martin Uecker <uecker@eecs.berkeley.edu>
7 #
8 # Compressed sensing parallel imaging reconstruction with temporal
9 # total-variation regularization for Siemens radial VIBE sequence
10 # with golden-angle sampling (GRASP).
11 #
12 set -e
13
14 # default settings
15 export SPOKES=21
16 export SKIP=0
17 export CALIB=400
18 export ITER=30
19 export REG=0.05
20 SCALE=0.6
21 LOGFILE=/dev/null
22 MAXPROC=4
23 MAXTHREADS=4
24
25 title=$(cat <<- EOF
26 (BART-)GRASP v0.2 (Berkeley Advanced Reconstruction Toolbox)
27 --- EXPERIMENTAL --- FOR RESEARCH USE ONLY ---
28 EOF
29 )
30
31 helpstr=$(cat <<- EOF
32 Compressed sensing parallel imaging reconstruction with temporal
33 total-variation regularization for Siemens radial VIBE sequence
34 with golden-angle sampling (GRASP).
35 This script requires the Berkeley Advanced Reconstruction Toolbox
36 version 0.2.08. Later versions may work.
37
38 -s spokes number of spokes per frame
39 -r lambda regularization parameter
40 -p maxproc max. number of slices processed in parallel
41 -t maxthreads max. number of threads per slice
42 -l logfile
43 -h help
44 EOF
45 )
46
47 usage="Usage: $0 [-h] [-s spokes] [-r lambda] <meas.dat> <output>"
48
49 echo "$title"
50 echo
51
52 while getopts "hl:s:p:t:r:" opt; do
53 case $opt in
54 s)
55 SPOKES=$OPTARG
56 ;;
57 r)
58 REG=$OPTARG
59 ;;
60 h)
61 echo "$usage"
62 echo
63 echo "$helpstr"
64 exit 0
65 ;;
66 l)
67 LOGFILE=$(readlink -f "$OPTARG")
68 ;;
69 p)
70 MAXPROC=$OPTARG
71 ;;
72 t)
73 MAXTHREADS=$OPTARG
74 ;;
75 \?)
76 echo "$usage" >&2
77 exit 1
78 ;;
79 esac
80 done
81
82 shift $((OPTIND - 1))
83
84
85 if [ $# -lt 2 ] ; then
86
87 echo "$usage" >&2
88 exit 1
89 fi
90
91
92 export PATH=$TOOLBOX_PATH:$PATH
93
94 input=$(readlink -f "$1")
95 output=$(readlink -f "$2")
96
97 if [ ! -e $input ] ; then
98 echo "Input file does not exist." >&2
99 echo "$usage" >&2
100 exit 1
101 fi
102
103 if [ ! -e $TOOLBOX_PATH/fft ] ; then
104 echo "\$TOOLBOX_PATH is not set correctly!" >&2
105 exit 1
106 fi
107
108
109 #WORKDIR=$(mktemp -d)
110 # Mac: http://unix.stackexchange.com/questions/30091/fix-or-alternative-for-mktemp-in-os-x
111 WORKDIR=`mktemp -d 2>/dev/null || mktemp -d -t 'mytmpdir'`
112 trap 'rm -rf "$WORKDIR"' EXIT
113 cd $WORKDIR
114
115 # start group for redirection of output to the logfile
116 {
117
118 # read TWIX file
119 twixread -A $input grasp
120
121 export READ=$(tail -n1 grasp.hdr | cut -f1 -d" ")
122 export COILS=$(tail -n1 grasp.hdr | cut -f4 -d" ")
123 export PHASES=$(($(tail -n1 grasp.hdr | cut -f2 -d" ") / $SPOKES))
124
125 export OMP_NUM_THREADS=$((MAXPROC * $MAXTHREADS))
126
127 # zero-pad
128 #flip $(bitmask 2) grasp grasp2
129 #resize 2 64 grasp2 grasp
130 #circshift 2 10 grasp grasp2
131 #fft -u $(bitmask 2) grasp2 grasp_hybrid
132 #rm grasp.* grasp2.*
133
134 # inverse FFT along 3rd dimension
135 fft -i -u $(bitmask 2) grasp grasp_hybrid
136 rm grasp.cfl grasp.hdr
137
138 SLICES=$(tail -n1 grasp_hybrid.hdr | cut -f3 -d" ")
139
140
141 # create trajectory with 400 spokes and 2x oversampling
142 traj -G -x$READ -y$CALIB r
143 scale $SCALE r rcalib
144
145 # create trajectory with 2064 spokes and 2x oversampling
146 traj -G -x$READ -y$(($SPOKES * $PHASES)) r
147 scale $SCALE r r2
148
149 # split off time dimension into index 10
150 reshape $(bitmask 2 10) $SPOKES $PHASES r2 rfull
151
152 # number of threads per slice
153 export OMP_NUM_THREADS=$MAXTHREADS
154
155 calib_slice()
156 {
157 # extract slice
158 slice 2 $1 grasp_hybrid grasp1-$1
159
160 # extract first $CALIB spokes
161 extract 1 $(($SKIP + 0)) $(($SKIP + $CALIB - 1)) grasp1-$1 grasp2-$1
162
163 # reshape dimensions
164 reshape $(bitmask 0 1 2 3) 1 $READ $CALIB $COILS grasp2-$1 grasp3-$1
165
166 # apply inverse nufft to first $CALIB spokes
167 nufft -i -t rcalib grasp3-$1 img-$1.coo
168 }
169
170 recon_slice()
171 {
172 # extract sensitivities for slice
173 slice 2 $1 sens sens-$1
174
175 # extract spokes and split-off time dim
176 extract 1 $(($SKIP + 0)) $(($SKIP + $SPOKES * $PHASES - 1)) grasp1-$1 grasp2-$1
177 reshape $(bitmask 1 2) $SPOKES $PHASES grasp2-$1 grasp1-$1
178
179 # move time dimensions to dim 10 and reshape
180 transpose 2 10 grasp1-$1 grasp2-$1
181 reshape $(bitmask 0 1 2) 1 $READ $SPOKES grasp2-$1 grasp1-$1
182 rm grasp2-$1.cfl grasp2-$1.hdr
183
184 # reconstruction with tv penality along dimension 10
185 pics -S -d5 -lv -u10. -r$REG -R$(bitmask 10) -i$ITER -t rfull grasp1-$1 sens-$1 i-$1.coo
186
187 # clean up temp files
188 rm *-$1.cfl *-$1.hdr
189 }
190
191 export -f calib_slice
192 export -f recon_slice
193
194 # loop over slices
195 seq -w 0 $(($SLICES - 1)) | xargs -i --max-procs=$MAXPROC bash -c "calib_slice {}"
196
197 # transform back to k-space and compute sensitivities
198 join 2 img-*.coo img
199 fft -u $(bitmask 0 1 2) img ksp
200
201 #ecalib -S -c0.8 -m1 -r20 ksp sens
202
203 # transpose because we already support off-center calibration region
204 # in dim 0 but here we might have it in 2
205 transpose 0 2 ksp ksp2
206 ecalib -S -c0.8 -m1 -r20 ksp2 sens2
207 transpose 0 2 sens2 sens
208
209 # loop over slices
210 seq -w 0 $(($SLICES - 1)) | xargs -i --max-procs=$MAXPROC bash -c "recon_slice {}"
211 #echo 20 | xargs -i --max-procs=$MAXPROC bash -c "recon_slice {}"
212
213 # join slices back together
214 join 2 i-*.coo $output
215
216 # generate dicoms
217 #for s in $(seq -w 0 $(($SLICES - 1))) ; do
218 # for p in $(seq -w 0 $(($PHASES - 1))) ; do
219 #
220 # slice 10 $p i-$s.coo i-$p-$s.coo
221 # toimg i-$p-$s.coo $output-$p$-$s.dcm
222 # done
223 #done
224
225 } > $LOGFILE
226
227 exit 0
228
229