%% Fourier phase retrieval with phase-uncertain mask
%% P random illuminations
%% Miao's oversampling with zero paddings
%% Wenjing Liao (wjliao@math.duke.edu)

clear all
close all
display('--------------------begin-------------------------')
%% parameters
i = sqrt(-1);
P = 2; % the number illuminations
Sample2 = 4; % total oversampling rate in two dimensions
Sample1 = sqrt(Sample2); % oversampling rate at each dimension
constmask = 1; % 1: use one uniform illumination
% 0: no uniform illumination, all random
brow = 1; % always set to 1
bcol = 1; % always set to 1
smoothmask = 1; % 0 HRM: High resolution mask
% 1 LRM: Low resolution mask
smoothobject= 0; % For complex-valued images
% if smoothobject == 1, images have low-resolution random
% phases
% if smoothobject == 0, images have high-resolution
% random phases
%% construct moliffier to generate LRM
ma = 4;
[xx , yy] = meshgrid(-ma+1:ma-1,-ma+1:ma-1);
Mollifier = exp(-ma^2./(ma^2-xx.^2)).*exp(-ma^2./(ma^2-yy.^2));
%%% normalize the moliffier
Mollifier = Mollifier/sum(sum(Mollifier));
%% object mollifier to generate images with low-resolution random phases when smoothmask = 1
mo = 5;
[xx , yy] = meshgrid(-mo+1:mo-1,-mo+1:mo-1);
OMollifier = exp(-mo^2./(mo^2-xx.^2)).*exp(-mo^2./(mo^2-yy.^2));
%%% normalize the moliffier
OMollifier = OMollifier/sum(sum(OMollifier));

%% measurement noise type and noise level
noisetype = 0; % 0 gaussian noise 1 poisson noise
noiselevel = 0; % external error

masknoise = 30; % this corresponds phase uncertainty in paper
% eg. masknoise=30 <=> \delta=0.3

testimage = 4; % 1 cameraman % 2 virus % 3 bull % 4 phantom %5 mandrill % 6 face
imagetype = -2; % 1 real positive % 2 complex positive
% -2 complex
ifdisplay = 1;

itermax = 200 + masknoise*10;
ifer = 200 + masknoise*10;
ifadaptive = 1; % adapt the termination of DRER and AER

%% Mask update rule
% maskupdate == 3
% standard rule for mask update in the paper ''Fourier phasing with phase-uncertain mask''
% if mask angle constraint is violated, do orthogonal projection to enforce the mask angle constraint
% maskupdate == 1
% alternative rule for mask update in the paper ''Fourier phasing with phase-uncertain mask''
% if mask angle constraint is violated, use the mask estimate
if imagetype == 1
maskupdate = 3;
elseif imagetype == 2
maskupdate = 1;
elseif imagetype == -2
maskupdate = 1;
end

%% loose support or tight support for images
supptype = 1; % 0 tight support 1 loose support
if supptype == 0
loose = 0;
elseif supptype == 1
if testimage == 1
loose = 13;
elseif testimage == 2
loose = 10;
elseif testimage == 3
loose = 10;
elseif testimage == 4
loose = 0;
elseif testimage == 5
loose = 0;
elseif testimage == 6
loose = 0;
end
end

%% load real positive images
if testimage == 1
Xsingle = imread('cameraman.tif');
elseif testimage == 2
Xsingle = imread('influenza.jpg');
Xsingle = rgb2gray(Xsingle);
elseif testimage == 3
Xsingle = imread('bull.jpg');
Xsingle = rgb2gray(Xsingle);
elseif testimage == 4
Xsingle = phantom('Modified Shepp-Logan',200);
elseif testimage == 5
Xsingle = imread('mandrill.jpg');
Xsingle = rgb2gray(Xsingle);
elseif testimage == 6
Xsingle = imread('face.png');
Xsingle = rgb2gray(Xsingle);
end
Xsingle = im2double(Xsingle); % image without zero paddings
[N1 , N2] = size(Xsingle);
if smoothobject == 0 % complex images with high-resolution random phases
if imagetype == 2
Xsingle = Xsingle.*(exp(i*unifrnd(0,pi/2,N1,N2)));
elseif imagetype == -2
Xsingle = Xsingle.*exp(i*unifrnd(0,pi*2,N1,N2));
end
elseif smoothobject == 1 % complex images with low-resolution random phases
if imagetype == 2
ImagePhase = 2*pi*rand(N1,N2);
TemPhase = exp(i*ImagePhase);
TemPhase = conv2(TemPhase,OMollifier,'same');
ImagePhase = 1/4*angle(TemPhase)+pi/4;
Xsingle = Xsingle.*exp(i*ImagePhase);
elseif imagetype == -2
ImagePhase = 2*pi*rand(N1,N2);
TemPhase = exp(i*ImagePhase);
TemPhase = conv2(TemPhase,OMollifier,'same');
ImagePhase = angle(TemPhase);
Xsingle = Xsingle.*exp(i*ImagePhase);
end
end

%% generate images with zero paddings to use the oversampling technique by
%% Miao, et. al.
N11 = round(Sample1 * (N1+2*loose));
N22 = round(Sample1 * (N2+2*loose));
Xdouble = zeros(N11,N22);
N1start = round(1 + (N11-N1)/2);
N1end = N1start + N1 - 1;
N2start = round(1 + (N22-N2)/2);
N2end = N2start + N2 - 1;

S1start = N1start - loose;
S1end = N1end + loose;
S2start = N2start - loose;
S2end = N2end + loose;

Xdouble(N1start:N1end,N2start:N2end) = Xsingle;

%% The following part generates exact mask: Diagonal
%% and mask estimate: Diagonal_noisy
% random entry at each block
if constmask == 1 % one uniform illumination
indexstart = 2;
elseif constmask == 0 % no uniform illumination, all random
indexstart = 1;
end
% the number of blocks
nbrow = ceil((S1end-S1start+1)/brow);
nbcol = ceil((S2end-S2start+1)/bcol);
% random phase
Diagonal = zeros(N11,N22,P);
Diagonal_noisy = zeros(N11,N22,P);
BlockPhase = zeros(nbrow,nbcol,P);
for k = indexstart : P
BlockPhase(:,:,k) = 2*pi*rand(nbrow,nbcol);
end
RandomPhase = zeros(S1end-S1start+1,S2end-S2start+1,P);
for k = 1 : nbrow
for l = 1 : nbcol
for t = 1 : P
RandomPhase((k-1)*brow+1:min(S1end-S1start+1,k*brow),(l-1)*bcol+1:min(S2end-S2start+1,l*bcol),t) = BlockPhase(k,l,t);
end
end
end
% if smoothmask == 1, create LRM
if smoothmask == 1
for k = indexstart : P
TemMask = exp(i*RandomPhase(:,:,k));
TemMask = conv2(TemMask,Mollifier,'same');
RandomPhase(:,:,k) = angle(TemMask);
end
end
Diagonal(S1start:S1end,S2start:S2end,:) = exp(i*RandomPhase);
% add noise to mask
PhaseError = zeros(S1end-S1start+1,S2end-S2start+1,P);
if masknoise > 0
for k = indexstart : P
if smoothmask == 0
PhaseError(:,:,k) = unifrnd(-pi*masknoise/100,pi*masknoise/100,S1end-S1start+1,S2end-S2start+1);
elseif smoothmask == 1
PhaseError(:,:,k) = 2*pi*rand(S1end-S1start+1,S2end-S2start+1);
TemMask = exp(i*PhaseError(:,:,k));
TemMask = conv2(TemMask,Mollifier,'same');
PhaseError(:,:,k) = masknoise/100*angle(TemMask);
end
end
RandomPhase_noisy = RandomPhase + PhaseError;
else
RandomPhase_noisy = RandomPhase;
end
Diagonal_noisy(S1start:S1end,S2start:S2end,:) = exp(i*RandomPhase_noisy);
%%
TEMD = Diagonal_noisy;
Diagonal_noisy = Diagonal;
Diagonal = TEMD;

%% noise-free high resolution random phase illumination in padding
for k = 1 : P
Diagonal(1:N11,1:S2start-1,k) = exp(i*2*pi*rand(N11,S2start-1));
Diagonal(1:N11,S2end+1:N22,k) = exp(i*2*pi*rand(N11,N22-S2end));
Diagonal(1:S1start-1,S2start:S2end,k) = exp(i*2*pi*rand(S1start-1,S2end-S2start+1));
Diagonal(S1end+1:N11,S2start:S2end,k) = exp(i*2*pi*rand(N11-S1end,S2end-S2start+1));
end
Diagonal_noisy(1:N11,1:S2start-1,:) = Diagonal(1:N11,1:S2start-1,:);
Diagonal_noisy(1:N11,S2end+1:N22,:) = Diagonal(1:N11,S2end+1:N22,:);
Diagonal_noisy(1:S1start-1,S2start:S2end,:) = Diagonal(1:S1start-1,S2start:S2end,:);
Diagonal_noisy(S1end+1:N11,S2start:S2end,:) = Diagonal(S1end+1:N11,S2start:S2end,:);

X = Xdouble;
%% Generate Fourier intensity measurements with exact mask: Diagonal
Y = zeros(N11,N22,P);
Y_abs_fro = 0;
for k = 1 : P
Y(:,:,k) = fft2(Diagonal(:,:,k).*X);
Y_abs_fro = Y_abs_fro + norm(Y(:,:,k),'fro')^2;
end
Y_abs_fro = sqrt(Y_abs_fro);
Y_abs = abs(Y);
Y_abs_noisy = zeros(size(Y_abs));

%% Add external noise to measurements
if noisetype == 0
sigma = noiselevel * Y_abs_fro / (100 * sqrt(2*P*N11*N22));
Enoise = zeros(N11,N22,P);
E_fro = 0;
for k = 1 : P
Enoise(:,:,k) = normrnd(0,sigma,N11,N22) + i * normrnd(0,sigma,N11,N22);
E_fro = E_fro + norm(Enoise(:,:,k),'fro')^2;
end
E_fro = sqrt(E_fro);
Y_noisy = Y + Enoise;
Y_abs_noisy = abs(Y_noisy);
elseif noisetype == 1
if noiselevel > 0
pscale = (100/noiselevel)^2*(sum(sum(sum(Y_abs))))/(Y_abs_fro^2);
for k = 1 : P
Enoise = poissrnd(pscale*Y_abs(:,:,k));
Y_abs_noisy(:,:,k) = Enoise / pscale;
end
E_fro = 0;
for k = 1 : P
E_fro = E_fro + norm(Y_abs(:,:,k)-Y_abs_noisy(:,:,k),'fro')^2;
end
E_fro = sqrt(E_fro);
elseif noiselevel == 0
Y_abs_noisy = Y_abs;
E_fro = 0;
end
end

%% starting point
Index = [N1 , N2 ; N11 , N22 ; S1start , S1end; S2start , S2end]; % parameters about images
X0 = generatestart(imagetype,Index);

%% Apply DRER followed by AER
Diagonal0 = Diagonal_noisy; % mask estimate
if masknoise == 0 % mask exactly known, apply DR followed by ER without mask update
[Xrec , IterRes , hiostep , erstep] = DR(Diagonal_noisy,Y_abs_noisy,imagetype,Index,itermax,ifer,X0,ifadaptive,ifdisplay);
else
[Xrec , Diagonalrec , IterRes, hiostep , erstep ] = DR_NoisyMask(Y_abs_noisy,imagetype,Index,itermax,ifer,X0,Diagonal0,masknoise,maskupdate,indexstart,ifadaptive,ifdisplay);
end

%% image error
theta = globalerror(X,Xrec);
spaceerr = norm(X-exp(i*theta)*Xrec,'fro')/norm(X,'fro');
fprintf('theta = %6.4f \n',theta)
fprintf('||X-eitheta*Xrec||/||X|| = %6.4f percent \n',100*spaceerr)

%% plot amplitudes of recovered images
figure
imshow(abs(Xrec(S1start:S1end,S2start:S2end)))
title(['image error = ' num2str(spaceerr*100) '%'])

%% Mask angle error
if masknoise > 0
Angle = angle(Diagonal(S1start:S1end,S2start:S2end,end));
Anglerec = angle(Diagonalrec(S1start:S1end,S2start:S2end,end));
AngleDiff = abs(Angle-Anglerec);
AngleDiff = min(AngleDiff,2*pi-AngleDiff);
figure
imshow(AngleDiff/pi)
title('Angle difference between exact and recovered mask')
end

%% Relative residual at each iteration
figure
plot(1:length(IterRes),IterRes,'.-')
axis([1 length(IterRes) 0 1])
title('relative residual at each iteration')
xlabel('iteration')
ylabel('relative residual')