-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathdetectorType.m
executable file
·301 lines (289 loc) · 12.2 KB
/
detectorType.m
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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
function [mtf, dqe, cf, readn, darkn] = detectorType(params)
%detectorType: Reads in the paratemres of already characterized detectors
% SYNOPSIS:
% [mtf, dqe, cf, readn, darkn] = detectorType(params)
%
% PARAMETERS:
% params: structure containing various input parameters (e.g. camera type, voltage, directory)
%
% OUTPUT:
% mtf: modulation transfer function
% dqe: detective quantum efficieny
% readn: readout noise
% darkn: dark current noise
% (C) Copyright 2013
% Quantitative Imaging Group Leiden University Medical Center
% Faculty of Applied Sciences Department of Molecular Cell Biology
% Delft University of Technology Section Electron Microscopy
% Lorentzweg 1 2300 RC Leiden
% 2628 CJ Delft
% The Netherlands
%
% Milos Vulovic
cam = params.cam.type;
Voltage = params.acquis.Voltage;
camdir = params.camdir;
s0= 'MTF is not available for this voltage: %03d keV. We will inter(extra)polate from the MTFs measured at %03d and %03d keV\n';
s01= 'MTF is not available for this voltage: %03d keV. We will extrapolate from the MTF measured at %03d keV\n';
s02= 'MTF is not available for this voltage: %03d keV. We will inter(extra)polate from the MTFs measured at %03d, %03d and %03d keV\n';
if (strcmp(cam,'Eagle4k') && Voltage/1000~=120)
s1=sprintf(s01, Voltage/1000, 120);
warning(s1);
elseif (strcmp(cam,'US4000') && (Voltage/1000~=120 && Voltage/1000~=200))
s1=sprintf(s0, Voltage/1000, 120, 200);
warning(s1);
elseif strcmp(cam,'Falcon') && (Voltage/1000~=200 && Voltage/1000~=300)
s1=sprintf(s0, Voltage/1000, 200, 300);
warning(s1);
elseif strcmp(cam,'US1000GIF') && (Voltage/1000~=80 && Voltage/1000~=200 && Voltage/1000~=300)
s1=sprintf(s02, Voltage/1000, 80, 200, 300);
warning(s1);
elseif strcmp(cam,'custom')
warning('This camera has not been characterized. Please provide a MTF and DQE. Or switch camera type to ''custom''');
end
if strcmp(cam,'ideal')
cf = 1000;
mtf = dip_image(ones(256,256));
dqe = mtf^2;
readn = 0;
darkn = 0;
elseif strcmp(cam,'perfect')
cf = 1000;
mtf = sinc(xx/size(xx,1))*sinc(yy/size(yy,1));
dqe = mtf^2;
readn = 0;
darkn = 0;
elseif strcmp(cam,'US1000GIF') % Gatan US1000. Data available for 80, 200 and 300 kV. For the other voltages interpolation is used
if Voltage == 80e3
cf = 28;
bz = load([camdir filesep 'MTF_US1000GIF_080']);
mtf = bz.mtf;
if exist([camdir filesep 'DQE_US1000GIF_080.mat'])
npsz = load([camdir filesep 'DQE_US1000GIF_080']);
dqe = npsz.dqe;
else
warning('DQE not available. Use only MTF.')
dqe = mtf^2;
end
readn = 3;
darkn = 0.11;
elseif Voltage == 200e3
cf = 15; %check
bz = load([camdir filesep 'MTF_US1000GIF_200']);
mtf = bz.mtf;
if exist([camdir filesep 'DQE_US1000GIF_200.mat'])
npsz = load([camdir filesep 'DQE_US1000GIF_200']);
dqe = npsz.dqe;
else
warning('DQE not available. Use only MTF.')
dqe = mtf^2;
end
readn = 3;
darkn = 0.11;
elseif Voltage == 300e3
cf = 6.6;
bz = load([camdir filesep 'MTF_US1000GIF_300']);
mtf = bz.mtf;
if exist([camdir filesep 'DQE_US1000GIF_300.mat'])
npsz = load([camdir filesep 'DQE_US1000GIF_300']);
dqe = npsz.dqe;
else
warning('DQE not available. Use only MTF.')
dqe = mtf^2;
end
readn = 3;
darkn = 0.11;
else
cf = (sqrt(80)*28*sqrt(200)*15*sqrt(300)*6.6)^(1/3)/sqrt(Voltage/1000);
bz = load([camdir filesep 'MTF_US1000GIF_080']);
mtf80 = bz.mtf;
bz = load([camdir filesep 'MTF_US1000GIF_200']);
mtf200 = bz.mtf;
bz = load([camdir filesep 'MTF_US1000GIF_300']);
mtf300 = bz.mtf;
mtf = (sqrt(80)*mtf80*sqrt(200)*mtf200*sqrt(300)*mtf300)^(1/3)/sqrt(Voltage/1000);
mtf = mtf/max(mtf);
if exist([camdir filesep 'DQE_US1000GIF_080.mat'])&& exist([camdir filesep 'DQE_US1000GIF_200.mat']) && exist([camdir filesep 'DQE_US1000GIF_300.mat'])
npsz = load([camdir filesep 'DQE_US1000GIF_080']);
dqe80 = npsz.dqe;
npsz = load([camdir filesep 'DQE_US1000GIF_200']);
dqe200= npsz.dqe;
npsz = load([camdir filesep 'DQE_US1000GIF_300']);
dqe300= npsz.dqe;
dqe = (sqrt(80)*dqe80*sqrt(200)*dqe200*sqrt(300)*dqe300)^(1/3)/sqrt(Voltage/1000);
else
warning('Noise power spectra not available. Use only MTF.')
dqe = mtf^2;
end
readn = 3;
darkn = 0.11;
end
elseif strcmp(cam,'Eagle4k')
if Voltage == 120e3
cf = 100; %check
bz = load([camdir filesep 'MTF_Eagle_120']);
mtf = bz.mtf;
if exist([camdir filesep 'DQE_Eagle_120.mat'])
npsz = load([camdir filesep 'DQE_Eagle_120']);
dqe = npsz.dqe;
else
warning('DQE not available. Use only MTF.')
dqe = mtf^2;
end
readn = 7;
darkn = 1.14;
else
cf = sqrt(120)*100/sqrt(Voltage/1000);
bz = load([camdir filesep 'MTF_Eagle_120']);
mtf120 = bz.mtf;
mtf = sqrt(120)*mtf120/sqrt(Voltage/1000);
mtf = mtf/max(mtf);
if exist([camdir filesep 'DQE_Eagle_120.mat'])
npsz = load([camdir filesep 'DQE_Eagle_120']);
dqe120= npsz.dqe;
dqe = sqrt(120)*dqe120/sqrt(Voltage/1000);
else
warning('DQE not available. Use only MTF.')
dqe = mtf^2;
end
readn = 7;
darkn = 1.14;
end
elseif strcmp(cam,'US4000')
if Voltage == 120e3
cf = 34;
bz = load([camdir filesep 'MTF_US4000_120']);
mtf = bz.mtf;
if exist([camdir filesep 'DQE_US4000_120.mat'])
npsz = load([camdir filesep 'DQE_US4000_120']);
dqe = npsz.dqe;
else
warning('DQE not available. Use only MTF.')
dqe = mtf^2;
end
readn = 3;
darkn = 0.11;
elseif Voltage == 200e3
cf = 23;
bz = load([camdir filesep 'MTF_US4000_200']);
mtf = bz.mtf;
if exist([camdir filesep 'DQE_US4000_200.mat'])
npsz = load([camdir filesep 'DQE_US4000_200']);
dqe = npsz.dqe;
else
warning('DQE not available. Use only MTF.')
dqe = mtf^2;
end
readn = 3;
darkn = 0.11;
else
cf = sqrt(sqrt(120)*34*sqrt(200)*23)/sqrt(Voltage/1000);
bz = load([camdir filesep 'MTF_US4000_120']);
mtf120 = bz.mtf;
bz = load([camdir filesep 'MTF_US4000_200']);
mtf200 = bz.mtf;
mtf = sqrt(sqrt(120)*mtf120*sqrt(200)*mtf200)/sqrt(Voltage/1000);
mtf = mtf/max(mtf);
if exist([camdir filesep 'DQE_US4000_120.mat'])&& exist([camdir filesep 'DQE_US4000_200.mat'])
npsz = load([camdir filesep 'DQE_US4000_120']);
dqe120 = npsz.dqe;
npsz = load([camdir filesep 'DQE_US4000_200']);
dqe200 = npsz.dqe;
dqe = sqrt(sqrt(120)*dqe120*sqrt(200)*dqe200)/sqrt(Voltage/1000);
else
warning('Noise power spectra not available. Use only MTF.')
dqe = mtf^2;
end
readn = 3;
darkn = 0.11;
end
elseif strcmp(cam,'FalconI')
if Voltage == 200e3
cf = 120;
bz = load([camdir filesep 'MTF_Falcon_200']);
mtf = bz.mtf;
npsz = load([camdir filesep 'DQE_Falcon_200']);
dqe = npsz.dqe;
readn = 3;
darkn = 0.11;
elseif Voltage == 300e3
cf = 100;
bz = load([camdir filesep 'MTF_Falcon_300']);
mtf = bz.mtf;
npsz = load([camdir filesep 'DQE_Falcon_300']);
dqe = npsz.dqe;
readn = 2;
darkn = 0.11;
else
cf = 100;
bz = load([camdir filesep 'MTF_Falcon_200']);
mtf200 = bz.mtf;
bz = load([camdir filesep 'MTF_Falcon_300']);
mtf300 = bz.mtf;
mtf = sqrt(sqrt(200)*mtf200*sqrt(300)*mtf300)/sqrt(Voltage/1000);
mtf = mtf/max(mtf);
if exist([camdir filesep 'DQE_Falcon_200.mat'])&& exist([camdir filesep 'DQE_Falcon_300.mat'])
npsz = load([camdir filesep 'DQE_Falcon_200']);
dqe200 = npsz.dqe;
npsz = load([camdir filesep 'DQE_Falcon_300']);
dqe300 = npsz.dqe;
dqe = sqrt(sqrt(200)*dqe200*sqrt(300)*dqe300)/sqrt(Voltage/1000);
else
warning('Noise power spectra not available. Use only MTF.')
dqe = mtf^2;
end
readn = 3;
darkn = 0.11;
end
elseif strcmp(cam,'K2')
if Voltage == 300e3
cf = 1;
bz = load([camdir filesep 'MTF_K2Summit_300']);
mtf = dip_image(bz.mtf);
npsz = load([camdir filesep 'DQE_K2Summit_300']);
dqe = dip_image(npsz.dqe);
readn = 0;
darkn = 0;
end
elseif strcmp(cam,'custom')
warning('Parameters of this camera are not known. Please characterize it with the supplementary software or simulate it as EMG (exponentially modified Gaussian) by setting the flag "params.cam.GenMTFasEMG" ');
if params.cam.GenMTFasEMG
mtf = SimMTFasEMG(params);
if params.cam.DQEflag
dqe = params.cam.dqe0*(mtf/mtf^params.cam.mtf_exp4nnps+eps)^2;
cf = 1000*sqrt(params.cam.dqe0./(params.acquis.dose_on_sample*(params.acquis.pixsize*1e10)^2));
else
dqe = 1+0*mtf;
cf = 1000*sqrt(1./(params.acquis.dose_on_sample*(params.acquis.pixsize*1e10)^2));
end
end
readn = params.cam.readn;
darkn = params.cam.darkn;
else
error('This camera type is not known\n')
end
% Avoid possible NaN
mtf(mtf==0) = 2;
minmtf=min(mtf);
mtf= mtf*(mtf~=2)+dip_image(minmtf*double((mtf==2)));
dqe(dqe==0) = 2;
mindqe=min(dqe);
dqe= dqe*(dqe~=2)+dip_image(mindqe*double((dqe==2)));
if isfield(params.cam, 'nnps_empirical') && ~strcmp(cam,'ideal') && ~strcmp(cam,'perfect')
if params.cam.nnps_empirical
dqe = params.cam.dqe0*(mtf/mtf^params.cam.mtf_exp4nnps+eps)^2;
else
dqe = params.cam.dqe0*(1+0*mtf);
end
if strcmp(cam,'custom')
cf = 1000*sqrt(params.cam.dqe0./(params.acquis.dose_on_sample*(params.acquis.pixsize*1e10)^2));
end
end
% MTF
mtf = resample(mtf,params.proc.N/256);
mtf1 = cut(mtf,params.proc.N/params.cam.bin);
mtf = resample(mtf1, params.cam.bin);
% DQE
dqe = resample(dqe,params.proc.N/256);
dqe1 = cut(dqe,params.proc.N/params.cam.bin);
dqe = resample(dqe1, params.cam.bin);