DictionaryBasedBackgroundEstimation/example.jl at main · IBIResearch/DictionaryBasedBackgroundEstimation · 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
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
##################################
# Load Packages and Activate Env #
##################################
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using MPIReco, Plots, LazyArtifacts

# Download data
datadir = artifact"data"
store = MDFDatasetStore(datadir)

##############
# Parameters #
##############

# Dataset store handling
study = getStudies(store, "BGDrift")[1]
exp = getExperiment(study, 1)

# Measurement data
filenameMeas = path(exp)
f = MPIFile(filenameMeas)

# System matrices
datadirSF = joinpath(calibdir(store), "1.mdf")
fSF = MPIFile(datadirSF)

# Background frames
bgFrames = 51:55

# Selected frames
frames = 131  # change this to e.g. 51:191 to reconstruct a range of frames

# Bounds for frequency selection
minFreq = 20e3 # change this to 40e3 to exclude the excitation frequency
maxFreq = 200e3

# Use only the first receive channel
recChannels = 1:1

# Solver parameters
λ = 1.0
iterations = 20

# Parameters for the joint estimation of BG and particle concentration
beta = (1/5)^(5)
dictSize = 10

# Put all parameters into a dictionary
params = Dict{Symbol,Any}(
  :frames => frames,
  :bgFrames => bgFrames,
  :bgCorrectionInternal => true,
  :λ => λ,
  :recChannels => recChannels,
  :iterations => iterations,
  :minFreq => minFreq,
  :maxFreq => maxFreq,
  :SNRThresh => -1,
  :enforceReal => true,
  :enforcePositive => true,
  :spectralLeakageCorrection =>false
 )

##########################
# Perform reconstruction #
##########################

@info "Reconstruction with static BG correction"
paramsStatic = deepcopy(params)
cStaticCorr = reconstruction(fSF,f; paramsStatic...).data.data

@info "Reconstruction with linear interpolated BG correction"
paramsInterp = deepcopy(params)
paramsInterp[:bgFramesPost] = 201:205
cDynCorr = reconstruction(fSF,f; paramsInterp...).data.data

@info "Reconstruction with joint estimation of BG and particle distribution"
paramsJoint = deepcopy(params)
paramsJoint[:bgDictSize] = dictSize
paramsJoint[:beta] = beta
cBGEst = reconstruction(fSF,f; paramsJoint...).data.data

###########################
## Visualize the results ##
###########################

# Windowing parameters
vmin = 0.0
vmax = 0.0006

p1 = heatmap(cStaticCorr[1,:,:,1,1], clim=(vmin,vmax), c=:viridis,
             title="Static", aspect_ratio=1 )
p2 = heatmap(cDynCorr[1,:,:,1,1], clim=(vmin,vmax), c=:viridis,
             title="Linear", aspect_ratio=1 )
p3 = heatmap(cBGEst[1,:,:,1,1], clim=(vmin,vmax), c=:viridis,
             title="Joint", aspect_ratio=1 )

plot(p1, p2, p3, layout = (1, 3), size=(1000,350))