OptimizedSystemMatrixSampling/smRecovery.jl at main · IBIResearch/OptimizedSystemMatrixSampling · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
using MPIReco, JLD, PyPlot

####################
# load system matrix
####################
# Dataset store handling
# uncomment the next two lines if running this script separately
# datadir = artifact"MDFStore"
# store = MDFDatasetStore(datadir)

# system matrix
datadirSF2 = joinpath(calibdir(store), "2.mdf")
fSF2 = MPIFile(datadirSF2)
nx,ny,nz = calibSize(fSF2)
S2 = reshape( getSystemMatrix(fSF2,bgCorrection=true), nx, ny, :)
snr2 = MPIReco.getSNRAllFrequencies(fSF2)

# select a slice and apply SNR thresholding
snr_idx2 = findall(x->x>3.0, vec(snr2))[2:end]
S2 = S2[:,:,snr_idx2]

########################
# undersampled data
########################
rf = 6
numSamp = div(nx*ny,6)
#poisson disk
idx_pd = load(datadir*"/samplingPatterns/pd_r$(rf).jld","idx")[1:numSamp]
# optimized patterns
# uncomment the next line if running this script without smSampling.jl
# idx_opt= load(datadir*"/samplingPatterns/samplingIdxOpt.jld","idx")[1:numSamp]

# perform undersampling
y_pd = zeros(ComplexF64, length(idx_pd),size(S2,3))
y_opt = zeros(ComplexF64, length(idx_opt),size(S2,3))
for k=1:size(S2,3)
  y_pd[:,k] .= vec(S2[:,:,k])[idx_pd]
  y_opt[:,k] .= vec(S2[:,:,k])[idx_opt]
end

############
# SMRecovery
############
params = Dict{Symbol,Any}()
params[:shape] = (nx,ny)      # size of reconstruction grid
params[:reg1] = "L1"          # type of regularization: L1
params[:sparseTrafo] = DCTOp(ComplexF64,(nx,ny),2)  # sparsifying transform
params[:ρ_l1] = 0.3
# params[:reg2] = "Nothing"
params[:iterationsInner] = 50 # maximu number of inner Split-Bregman iterations
params[:iterations] = 10      # number of outer split Bregman iterations
params[:relTol] = 1.e-3       # relative stopping tolerance for inner Split Bregman iteration
params[:absTol] = 1.e-4       # absolute stopping tolerance for inner Split Bregman iteration

params[:λ_l1] = 0.4^6         # ℓ1-regularization parameter
S2_pd = reshape( smRecovery(y_pd,idx_pd,params), nx, ny, :)

params[:λ_l1] = 0.4^5         # ℓ1-regularization parameter
S2_opt = reshape( smRecovery(y_opt,idx_opt,params), nx, ny, :)

################################
# plot some frequency components
################################
freqs=[24,35,38]
figure("recovered SMs", figsize=(6,6))
subplot(3,3,1); imshow(abs.(S2[:,:,freqs[1]])); title("k=$(snr_idx2[freqs[1]])"); ylabel("ref"); xticks([]); yticks([])
subplot(3,3,2); imshow(abs.(S2[:,:,freqs[2]])); title("k=$(snr_idx2[freqs[2]])"); xticks([]); yticks([])
subplot(3,3,3); imshow(abs.(S2[:,:,freqs[3]])); title("k=$(snr_idx2[freqs[3]])"); xticks([]); yticks([])
subplot(3,3,4); imshow(abs.(S2_pd[:,:,freqs[1]])); ylabel("PD"); xticks([]); yticks([])
subplot(3,3,5); imshow(abs.(S2_pd[:,:,freqs[2]])); xticks([]); yticks([])
subplot(3,3,6); imshow(abs.(S2_pd[:,:,freqs[3]])); xticks([]); yticks([])
subplot(3,3,7); imshow(abs.(S2_opt[:,:,freqs[1]])); ylabel("Opt"); xticks([]); yticks([])
subplot(3,3,8); imshow(abs.(S2_opt[:,:,freqs[2]])); xticks([]); yticks([])
subplot(3,3,9); imshow(abs.(S2_opt[:,:,freqs[3]])); xticks([]); yticks([])
tight_layout()