function [Xrec , iterres , iterchange] = HIOER(Diagonal,Y_abs_noisy,imagetype,Index,beta,X0,itermax,ifer)

%% parameters
i            = sqrt(-1);
N11          = Index(2,1);
N22          = Index(2,2);
N1start   = Index(3,1);
N1end     = Index(3,2);
N2start   = Index(4,1);
N2end     = Index(4,2);
P         = size(Diagonal,3);
N1        = N1end - N1start + 1;
N2        = N2end - N2start + 1;

ifadaptivehio = 1;
increase = 0;


%% initiation
Xk           = X0;      % Xk
XkP          = X0;      % Pf{Xk}
iterres      = [];      % record residual at each iteration
iterchange   = [];      % record ||XkP{k}-XkP{k-1}||/||XkP{k-1}||

Y_abs_noisy_fro = 0;    % Eucledean norm of Y_abs_noisy
for k = 1 : P
    Y_abs_noisy_fro           = Y_abs_noisy_fro + norm(Y_abs_noisy(:,:,k),'fro')^2;
end
Y_abs_noisy_fro = sqrt(Y_abs_noisy_fro);


%% HIO iteration
for iter = 1 : itermax  
    TemMatrix = XkP;     % record XkP from the previous iteration
    %% compute Pf{Xk}
    XkP   = Xk;
    for k = 1 : P
        FXk  = fft2(Diagonal(:,:,k).*XkP);
        FXkP = Y_abs_noisy(:,:,k).*exp(i*angle(FXk)); %*(FXk./abs(FXk));
        XkP  = ifft2(FXkP)./Diagonal(:,:,k);
    end
    %% apply HIO step   
    new     = zeros(N11,N22);
    %% apply HIO on zero paddings
    new(:,1:N2start-1)              = Xk(:,1:N2start-1) - beta * XkP(:,1:N2start-1);
    new(:,N2end+1:N22)              = Xk(:,N2end+1:N22) - beta * XkP(:,N2end+1:N22);
    new(1:N1start-1,N2start:N2end)  = Xk(1:N1start-1,N2start:N2end) - beta * XkP(1:N1start-1,N2start:N2end);
    new(N1end+1:N11,N2start:N2end)  = Xk(N1end+1:N11,N2start:N2end) - beta * XkP(N1end+1:N11,N2start:N2end);
    %% apply HIO in support
    if imagetype == 1
        newimag = imag(Xk(N1start:N1end,N2start:N2end)) - beta * imag(XkP(N1start:N1end,N2start:N2end));
        newreal = real(XkP(N1start:N1end,N2start:N2end));
        for m = 1 : N1
           for l = 1 : N2
               if newreal(m,l) < 0
                   newreal(m,l) = real(Xk(m+N1start-1,l+N2start-1)) - beta * real(XkP(m+N1start-1,l+N2start-1));
               end
           end
        end
        new(N1start:N1end,N2start:N2end) = newreal + i * newimag; 
        Xk             = new;
    elseif imagetype == 2
        newreal = real(XkP(N1start:N1end,N2start:N2end));
        newimag = imag(XkP(N1start:N1end,N2start:N2end));
        for m = 1 : N1
           for l = 1 : N2
               if newreal(m,l) < 0
                   newreal(m,l) = real(Xk(m+N1start-1,l+N2start-1)) - beta * real(XkP(m+N1start-1,l+N2start-1));
               end
               if newimag(m,l) < 0
                   newimag(m,l) = imag(Xk(m+N1start-1,l+N2start-1)) - beta * imag(XkP(m+N1start-1,l+N2start-1));
               end
           end
        end
        new(N1start:N1end,N2start:N2end) = newreal + i * newimag; 
        Xk             = new;
    elseif imagetype == -2
        new(N1start:N1end,N2start:N2end) = XkP(N1start:N1end,N2start:N2end); 
        Xk             = new;
    end
    %% record residual
    Xtem = zeros(N11,N22);
    Xtem(N1start:N1end,N2start:N2end) = projection(XkP(N1start:N1end,N2start:N2end),imagetype);
    intensityerr = 0;
    for k = 1 : P
        intensityerr = intensityerr + norm(abs(fft2(Diagonal(:,:,k).*Xtem))-Y_abs_noisy(:,:,k),'fro')^2;
    end
    intensityerr = sqrt(intensityerr) / Y_abs_noisy_fro;
    iterres(iter)  = intensityerr;
    %% record the change of XkP
    iterchange(iter) = norm(TemMatrix-XkP,'fro')/norm(TemMatrix,'fro');
    %% display data
    fprintf('%6.0f residual = %6.4f p change = %6.6f p \n',iter,100*intensityerr,100*iterchange(iter))
    %% Apply an adaptive rule for the termination of HIO
    %% The adaptive rule should vary from noiseless to noisy case
    %% The following rule is suitable for real and complex positive images in noisy case 
    %% For complex images, it is better not to terminate HIO early.
    if ifadaptivehio == 1
        if (imagetype==1) || (imagetype==2)
            %% stop if the residual doesn't drop in 5 steps
            if iterres(end) >= iterres(max(end-1,1))
                increase = increase + 1;
            else
                increase = 0;
            end
            if increase >= 5 && iter >= 10
                break;
            end
        end
    end
    %% Another adaptive rule which is suitable for noiseless case
% %     if iterchange(iter) < 0.01
% %         break
% %     end
end

hiostep  = length(iterres);
Xk       = XkP;   % the output of HIO is XkP


%% Apply ER
for iter = 1 : ifer
    TemMatrix = Xk;  % record Xk from the previous iteration 
    %% compute Pf{Xk}
    XkP   = Xk;
    for k = 1 : P
        FXk  = fft2(Diagonal(:,:,k).*XkP);
        FXkP = Y_abs_noisy(:,:,k).*(FXk./abs(FXk));
        XkP  = ifft2(FXkP)./Diagonal(:,:,k);
    end
    %% Po Pf{Xk}
    Xk  = zeros(N11,N22);
    Xk(N1start:N1end,N2start:N2end) = projection(XkP(N1start:N1end,N2start:N2end),imagetype);
    %% record residual
    intensityerr = 0;
    for k = 1 : P
        intensityerr = intensityerr + norm(abs(fft2(Diagonal(:,:,k).*Xk))-Y_abs_noisy(:,:,k),'fro')^2;
    end
    intensityerr = sqrt(intensityerr) / Y_abs_noisy_fro;
    iterres(hiostep+iter) = intensityerr;
    %% record the change of Xk
    iterchange(hiostep+iter) = norm(TemMatrix-Xk,'fro')/norm(TemMatrix,'fro');
    %% display data 
    fprintf('%6.0f residual = %6.4f p change = %6.6f p \n',iter,100*intensityerr,100*iterchange(hiostep+iter))
    %% Adaptive rule for ER in noiseless case
    %% stop if ||XkP{k}-XkP{k-1}||/||XkP{k-1}|| < 0.01%
%     if iterchange(hiostep+iter) < 0.0001
%         break;
%     end
end
     
    
Xrec   = Xk;
