forked from Waller-Lab/NOVOCGH
-
Notifications
You must be signed in to change notification settings - Fork 0
/
function_FunObj_dual.m
39 lines (36 loc) · 1.24 KB
/
function_FunObj_dual.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
function [loss, df ] = function_FunObj_dual( phase, source, z, Nx, Ny, thresholdh, thresholdl, maskFun, kickmaskFun,fresnelKernelFun, useGPU)
% Computes loss and gradient.
% This function is called by Matlab's fmincon library.
if useGPU
df = zeros(Nx, Ny, 'gpuArray');
phase = gpuArray(phase);
source = gpuArray(source);
else
df = zeros(Nx, Ny);
end
loss = 0;
phase = reshape(phase, [Nx, Ny]);
objectField = source.*exp(1i * phase);
for i = 1 : numel(z)
HStack = fresnelKernelFun(i,i);
mask = maskFun(i,i);
kickmask = kickmaskFun(i,i);
imagez = fftshift(fft2(objectField .* HStack));
imageInten = abs(imagez.^2);
maskh = mask .* (imageInten < thresholdh);
maskl = kickmask .* (imageInten > thresholdl);
diffh = maskh .* (imageInten - thresholdh);
diffl = maskl .* (imageInten - thresholdl);
temph = imagez.*diffh;
temph = conj(HStack).*(Nx*Ny*ifft2(ifftshift(temph)));
templ = imagez.*diffl;
templ = conj(HStack).*(Nx*Ny*ifft2(ifftshift(templ)));
%Compute losses
loss = loss + sum(sum(diffh.^2 + diffl.^2));
%Compute gradient
df = df + temph + templ;
end
dfphase = source.*(- real(df).*sin(phase) + imag(df) .* cos(phase));
df = gather(real(dfphase(:)));
loss = gather(real(loss));
end