-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfig_wiener.m
executable file
·121 lines (101 loc) · 4.34 KB
/
fig_wiener.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
% Draws graphs for the illustration of deconvolution
clear all
addpath('helper');
addpath('reconstruction');
N=1024;
objwidth=80;
objoffset=-150;
objdip=0.5;
objsum=1e6;
refsum=5e5;
refwidth=43;
refoffset=-objoffset;
noiselvl=1e-3;
normalize=@(in,ref)in./sum(in).*sum(ref);
ft =@(in) fftshift( fft(ifftshift(in)));
ift=@(in) fftshift(ifft(ifftshift(in)));
gaussFilter = @(in,sigma) normalize(ift(ft(exp( (-(N/2-(1:N)+1).^2) / (2 * sigma ^ 2))).*ft(in)),in);
ramp=linspace(1,objdip,floor(objwidth/4)+1);
rampoffset=objwidth-floor(objwidth/4)*4;
ramps=[ramp(1:end-1),fliplr(ramp(2:end)),ones(1,rampoffset),ramp(2:end-1),fliplr(ramp(1:end))];
obj=ramps;
obj=objsum*obj./sum(obj);
objpad=[zeros(1,ceil((N-objwidth)/2)+objoffset),obj,zeros(1,floor((N-objwidth)/2)-objoffset)];
objpad=gaussFilter(objpad,2);
ref=ones(1,refwidth);
ref=refsum*ref./sum(ref);
refpad=[zeros(1,ceil((N-refwidth)/2)+refoffset),ref,zeros(1,floor((N-refwidth)/2)-refoffset)];
refCentered=pad2size(ref,size(objpad));
refCentered=gaussFilter(refCentered,3);
refpad=gaussFilter(refpad,3);
signal=abs(objpad+refpad);
scatter=(abs(ft(signal)).^2);
noise=@(rel) (noiselvl*rel*(2*rand(1,N)-1));
holo=(ift(scatter));
n=noise(max(holo));
holo=holo+n;
cutsize=2*objwidth;
auto=pad2size(holo(1+end/2-objwidth:end/2+objwidth),[1,N]);
cross=pad2size(holo(1+end/2-abs(refoffset-objoffset)-cutsize:end/2-abs(refoffset-objoffset)+cutsize),[1,N]);
refcut=pad2size(refCentered(1+end/2-cutsize:end/2+cutsize),[1,N]);
objcut=pad2size(signal(1+end/2+objoffset-objwidth:1+end/2+objoffset+objwidth),size(cross));
autoref=ift(abs(ft(refcut)).^2);
deconv=abs(wiener(cross,refcut,abs(ft(n)).^2));
deconv1=abs(wiener(cross,refcut,abs(ft(n)).^2,ft(auto)));
deconv2=deconvwnr(cross,refcut,(ift(abs(ft(n)).^2)),auto);
deconv3=deconvwnr(cross,refcut,(ift(abs(ft(n)).^2)),ift(abs(ft(objcut)).^2));
directdeconv=(ift(ft(cross)./(1e-5+ft(refcut))));
mse=@(in)mean((objcut(:)-in(:)).^2);
deconvt=@(w)deconvwnr(cross,refcut,w);
options = optimset('TolFun',1e-18,'TolX',1e-8);
o=fminsearch(@(w)mse(deconvt(w)),1,options);
deconv4=deconvt(o);
figure(1);clf
blue=[48,38,131]./255;
green=[0,150,64]./255;
red=[227,5,19]./255;
yellow=[249,178,51]./255;
subplot(321);hold on;axis off;
plot(abs(objpad),'Color',green);
plot(abs(refpad),'Color',blue);
legend({'Objekt','Referenz'},'Interpreter','latex');
legend boxoff
zeroline = refline([0 0]);zeroline.Color = 'black';
subplot(322);hold on;axis off;
plot(log10(abs(ft(objcut))),'Color',green);
plot(log10(abs(ft(refcut))),'Color',blue);
xlim([N/8,7*N/8]);
a=ylim;ylim([a(1),a(2)*1.5]);
legend({'$\log \left|\mathcal{F}\left(\textrm{Objekt}\right)\right|$','$\log \left|\mathcal{F}\left(\textrm{Referenz}\right)\right|$'},'Interpreter','latex');
legend boxoff
zeroline = refline([0 0]);zeroline.Color = 'black';
subplot(323);hold on;axis off;
plot(abs(holo),'Color',red);
legend({'$\mathcal{F}^{-1}\left({\left|\mathcal{F}\left(\textrm{Objekt}+\textrm{Referenz}\right)\right|}^2\right)$'},'Interpreter','latex');
legend boxoff
zeroline = refline([0 0]);zeroline.Color = 'black';
subplot(324);hold on;axis off;
plot(log10(abs(ft(cross))),'Color',red);
plot(log10(abs(ft(refcut))),'Color',blue);
plot(log10(abs(ft(objcut))),'Color',green);
legend({'$\log \left|\mathcal{F}\left(\textrm{Kreuzkorrelation}\right)\right|$','$\log \left|\mathcal{F}\left(\textrm{Referenz}\right)\right|$','$\log \left|\mathcal{F}\left(\textrm{Objekt}\right)\right|$'},'Interpreter','latex');
legend boxoff
xlim([N/8,7*N/8]);
a=ylim;ylim([a(1),a(2)*1.5]);
zeroline = refline([0 0]);zeroline.Color = 'black';
subplot(325);axis off; hold on;
plot(abs(directdeconv./max(directdeconv).*max(objcut)),'Color',red);
plot(abs(deconv4),'Color',yellow);
legend({'direkte Entfaltung','Wiener Entfaltung'},'Interpreter','latex');
legend boxoff
zeroline = refline([0 0]);zeroline.Color = 'black';
subplot(326);hold on;axis off;
plot(log10(abs(ft(objcut))),'Color',green);
plot(log10(abs(ft(cross)))-log10(abs(ft(refcut))),'Color',red)
plot(log10(abs(ft(deconv4))),'Color',yellow);
legend({'$\log \left|\mathcal{F}\left(\textrm{Objekt}\right)\right|$','$\log \left|\mathcal{F}\left(\textrm{direkte Entfaltung}\right)\right|$','$\log \left|\mathcal{F}\left(\textrm{Wiener Entfaltung}\right)\right|$'},'Interpreter','latex');
legend boxoff
xlim([N/8,7*N/8]);
a=ylim;ylim([a(1),a(2)*1.5]);
zeroline = refline([0 0]);zeroline.Color = 'black';
print('Tex\images\src\wiener.pdf','-dpdf','-r600');