-
Notifications
You must be signed in to change notification settings - Fork 6
/
function_globalGS.m
54 lines (51 loc) · 1.49 KB
/
function_globalGS.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
%%% This is the implemetnation of GLOBAL Gerchberg Saxton algorithm.
%%% Optional, System.GSoffset hyperparameter
function [ GS ] = function_globalGS(System, HStacks, Masks )
if System.verbose == 1
disp('Global Gerchberg-Saxton hologram computation begins...');
tic;
end;
ims = Masks;
[NX,NY,NZ] = size(Masks);
if System.useGPU == 1
imageInten = zeros(NX,NY,NZ, 'gpuArray');
else
imageInten = zeros(NX,NY,NZ);
end
vv = System.verbose;
System.verbose = 0;
[ Superposition ] = function_Superposition( System,HStacks,Masks);
System.verbose = vv;
if System.verbose == 1
fprintf('Completed cycles #');
end;
Superposition.phase = mod(Superposition.phase, 2*pi) - pi;
im = System.source.*exp(1i * Superposition.phase);
for n = 1:System.maxiter
if System.verbose == 1
if mod(n, 10) == 0
fprintf(int2str(n));fprintf(',')
end
if n == System.maxiter;
fprintf(int2str(n));
end
end
index = 1:NZ;
tempim = 0 * im;
for i = index
imagez = fftshift(fft2(im .* HStacks(:,:,i)));
target = ims(:,:,i) + System.GSoffset;
imagez = sqrt(target) .* exp(1i * angle(imagez));
tempim = tempim + ifft2(ifftshift(imagez))./HStacks(:,:,i);
end
im = System.source.*exp(1i * angle(tempim));
end
%phase = gather(angle(im));
if System.verbose == 1
t = toc;
disp(' ');
disp(['Global Gerchberg-Saxton - Completed in ' int2str(t) ' seconds !']);
end;
GS.hologram = im;
GS.phase = gather(angle(GS.hologram));
end