reconstruction_refractive.m 2.25 KB
Newer Older
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
function reco = reconstruction_refractive(p)
% suppguess add the regions above and below the wire to the reconstruction.
% This is again for mitigation of boundary artifacts.
load('supp_guess.mat');
%% define the  projectors for the algorithm
% Plain initialisation for the guess with zeros in imaginary and real part.
% This is equivalent to vaccum.
guess = gpuArray(complex(zeros(p.recSize.*p.scale_fac, 'single')));
% projectors in the sample plane: refractive constraints and self refining
% support
PRef = PRefractive(p.refractiveSet);
PSuppAdapt = PSupportAdaption(gpuArray(zeros(p.recSize)), p.suppAdaptSet);
% projectors can be concatuated
PS = PRef * PSuppAdapt;

% magnitude constraint
p.PMset.applicationMask = [];
PM = PFresnelMagnitude(gpuArray(single(sqrt(p.holo))), p.F_calc, p.PMset);

%% 1st pass - find a support
starttime = datetime('now');
p.algSet.iterations = 1000;
projAlg = projectionAlgorithm(PM, PS, guess, p.algSet);
reco = projAlg.execute(p.algSet.iterations);

%% we refinie the support from the final reconstruchtion
seq_conv_hist{1} = reco.stats.convergenceHistory;
if(1)
p.supp = real(reco.x) < -0.2;
p.supp = padToSize(cropToCenter(p.supp, p.data_size), p.recSize, 0);
se = strel('disk', 70);
p.supp = p.supp | imdilate(supp, se);
p.supp = bwareafilt(p.supp, 1, 'largest');

se = strel('disk', 80);
p.supp = imerode(p.supp, se);
se = strel('disk', 80);
p.supp = imdilate(p.supp, se);

p.supp = imfill(p.supp, 'holes');
p.supp = bwareafilt(p.supp, 1, 'largest');
p.supp = imgaussfilt(single(p.supp), 10/2.35);
figure(1002)
imagesc(p.supp)
end
%% 2nd pass - keep support fixed. just iterate
if(1)
% define support Projector with constant support
PSupp = PSupport(gpuArray(p.supp), 0);
% set refractive constraints
p.refractiveSet.minReal = -80;
p.refractiveSet.minImag = -0.5e-3;
p.refractiveSet.maxImag = 500e-3;
% parameters for nesterov 
p.algSet.eta = 1.2;
p.algSet.gamma = 0.99;
PRef = PRefractive(p.refractiveSet);

PS = PRef * PSupp;% * PStat;
PM = PFresnelMagnitude(gpuArray(single(sqrt(p.holo))), p.F_calc, p.PMset);

p.algSet.iterations = 10000;
p.algSet.addMomentum = 1;
projAlg = projectionAlgorithm(PM, PS, guess, p.algSet);
reco = projAlg.execute(p.algSet.iterations);
seq_conv_hist{2} = reco.stats.convergenceHistory;
end

%% done
starttime - datetime('now')
end