@@ -27,7 +27,7 @@ please cite that paper.
27
27
# Packages needed here.
28
28
29
29
# # using Unitful: s
30
- using Plots; default (markerstrokecolor= :auto , label= " " )
30
+ using Plots; cgrad, default (markerstrokecolor= :auto , label= " " )
31
31
using MIRT: Afft, Asense, embed
32
32
using MIRT: pogm_restart, poweriter
33
33
using MIRTjim: jim, prompt
@@ -90,15 +90,16 @@ based on the
90
90
91
91
# ## Read data
92
92
if ! @isdefined (data)
93
- url = " https://web.eecs.umich.edu/~fessler/irt/reproduce/19/lin-19-edp/data/"
93
+ # src url = "https://web.eecs.umich.edu/~fessler/irt/reproduce/19/lin-19-edp/data/"
94
+ url = " https://github.com/JeffFessler/MIRTdata/raw/main/mri/lin-19-edp/"
94
95
dataurl = url * " cardiac_perf_R8.mat"
95
96
data = matread (Downloads. download (dataurl))
96
97
xinfurl = url * " Xinf.mat"
97
98
Xinf = matread (Downloads. download (xinfurl))[" Xinf" ][" perf" ] # (128,128,40)
98
99
end ;
99
100
100
101
# Show converged image as a preview:
101
- jim (Xinf, L "\m athrm{Converged\ image\ } X_∞" )
102
+ pinf = jim (Xinf, L "\m athrm{Converged\ image\ sequence } X_∞" )
102
103
103
104
# Organize k-space data:
104
105
if ! @isdefined (ydata0)
@@ -117,9 +118,22 @@ if !@isdefined(samp)
117
118
end
118
119
kx = - (nx÷ 2 ): (nx÷ 2 - 1 )
119
120
ky = - (ny÷ 2 ): (ny÷ 2 - 1 )
120
- jim (kx, ky, samp, " Sampling patterns for $nt frames" ; xlabel= L " k_x" , ylabel= L " k_y" )
121
+ psamp = jim (kx, ky, samp, " Sampling patterns for $nt frames" ;
122
+ xlabel= L " k_x" , ylabel= L " k_y" )
121
123
end
122
124
125
+ #=
126
+ Are all k-space rows are sampled in one of the 40 frames?
127
+ Sadly no.
128
+ The 10 blue rows shown below are never sampled.
129
+ A better sampling pattern design
130
+ could have avoided this issue.
131
+ =#
132
+ samp_sum = sum (samp, dims= 3 )
133
+ color = cgrad ([:blue , :black , :white ], [0 , 1 / 2 nt, 1 ])
134
+ pssum = jim (kx, ky, samp_sum; xlabel= " kx" , ylabel= " ky" ,
135
+ color, clim= (0 ,nt), title= " Number of sampled frames out of $nt " )
136
+
123
137
# Prepare coil sensitivity maps
124
138
if ! @isdefined (smaps)
125
139
smaps_raw = data[" b1" ] # raw coil sensitivity maps
@@ -130,7 +144,7 @@ if !@isdefined(smaps)
130
144
smaps = smaps_raw ./ ssos_raw
131
145
ssos = ssos_fun (smaps)
132
146
@assert all (≈ (1 ), ssos)
133
- jim (smaps, " Normalized |coil maps| for $nc coils" )
147
+ pmap = jim (smaps, " Normalized |coil maps| for $nc coils" )
134
148
end
135
149
136
150
@@ -156,7 +170,7 @@ because Xinf was reconstructed
156
170
using this regularizer!
157
171
=#
158
172
tmp = TF * Xinf
159
- jim (tmp, " |Temporal FFT of Xinf|" )
173
+ ptfft = jim (tmp, " |Temporal FFT of Xinf|" )
160
174
161
175
162
176
#=
211
225
212
226
# Check scale factor of Xinf. (It should be ≈1.)
213
227
tmp = A * Xinf
214
- scale = dot (tmp, ydata) / norm (tmp)^ 2 # 1.009 ≈ 1
228
+ scale0 = dot (tmp, ydata) / norm (tmp)^ 2 # 1.009 ≈ 1
215
229
216
230
# Crude initial image
217
231
L0 = A' * ydata # adjoint (zero-filled)
@@ -221,7 +235,7 @@ L0 = A' * ydata # adjoint (zero-filled)
221
235
S0 = zeros (ComplexF32, nx, ny, nt)
222
236
X0 = cat (L0, S0, dims= ndims (L0)+ 1 ) # (nx, ny, nt, 2) = (128, 128, 40, 2)
223
237
M0 = AII * X0 # L0 + S0
224
- jim (M0, " |Initial L+S via zero-filled recon|" )
238
+ pm0 = jim (M0, " |Initial L+S via zero-filled recon|" )
225
239
226
240
227
241
#=
0 commit comments