%% import data
addpath src;
addpath data;
load('data\cxidb-28_result3.mat');
[N1,N2] = size(dp);
dp(stopper)=-1;

%% GPS algorithm
clear input;

input.isreal      = false;
input.diffpat     = dp;
input.initial     = rand(N1,N2).*support;%u_initial;
input.support     = support;
input.sigma      = 0.01; 
input.iterations  = 4000;
[u2D, u2D_best] = GPS_algorithm(input);

%% HIO algorithm with averaging
u = rand(N1,N2).*support;
u_best  = u;
u_best0 = u;
err_best2 = 1;
err_best  = 1;
u_sum=0;
non_stopper = ~stopper;
no_supp = ~support;
norm_dp1 = norm ( dp(~stopper), 1 );
norm_dp2 = norm ( dp(~stopper), 2 );

iterations=4000;

for t=1:iterations
    Fu = fftshift(fft2(ifftshift(u)));
    diff = abs(Fu(non_stopper)) - dp(non_stopper);
    z_stopper = Fu(stopper);
    Fu = dp .* exp(1i*angle(Fu));
    Fu(stopper) = z_stopper;

    u_new = fftshift(ifft2(ifftshift(Fu)));


    %u_new(no_supp) = 0;
    u_new(no_supp) = u(no_supp) - 0.9*u_new(no_supp);
    %u_new = real(u_new);
    u = u_new;


    %error
    %
    Fu2  = fftshift(fft2(ifftshift(u.*support)));
    diff_z  = abs(Fu2(non_stopper))  - dp(non_stopper);
    err2 = norm(diff_z, 'fro') / norm_dp2;
    if err2<err_best2
        err_best2 =err2;
        u_best=u;
    end
    %}

    if t>3800
        u_sum=u_sum+u;
    end


    if mod(t,10)==0
        figure(4);
        rec = u_new.*support;
        img( rec, ['mag ',num2str(t)], 'colormap','gray')
        drawnow
        err1 = norm(diff_z(:),1) / norm_dp1;
        fprintf('%d. err1=%.3f, err2=%.3f\n',t,err1, err2);
    end
end
u_avg = u_sum/200;

%% result
u_avg_c = u_avg(360 +(1:128), 384+(1:256))';
u2D_c = u2D_best(360 +(1:128), 384+(1:256))';
u_c = u_best(360 +(1:128), 384+(1:256))';
figure(1); img(u_c,'HIO',u_avg_c,'HIO+avg', u2D_c,'GPS','colormap','gray', 'caxis',[0,7])

%%

%figure(2); img( u2D_c,'GPS','colormap','gray', 'caxis',[0,7])

%figure(3); img(u_c,'shrinkwrap', u2D_c,'GPS','colormap','gray', 'caxis',[0,7])





