function [Xrec , Diagonalk , iterres , hiostep , erstep ] = DR_NoisyMask(Y_abs_noisy,imagetype,Index,itermax,ifer,X0,Diagonal0,masknoise,maskupdate,indexstart,ifadaptive,ifdisplay)

%% 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(Diagonal0,3);

anglerange = masknoise/100*pi;
delta      = masknoise/100;

Angle0     = angle(Diagonal0);
%% initiation
Xk           = X0;
Imagek       = zeros(N11,N22);
Imagek(N1start:N1end,N2start:N2end) = projection(Xk(N1start:N1end,N2start:N2end),imagetype);
Y_abs_fro    = 0;
for k = 1 : P
    Y_abs_fro = Y_abs_fro + norm(Y_abs_noisy(:,:,k),'fro')^2;
end
Y_abs_fro    = sqrt(Y_abs_fro);

Diagonalk    = Diagonal0;

%% residual
iterres      = [];
%% for adaptive DRER
increase = 0;

%% DRER
for iter = 1 : itermax
    %% Pre
    Xpre = Xk;
    %% Po
    XkO = Imagek;
    %% Pf
    XkF = Xk;
    for k = 1 : P
        FXk  = fft2(Diagonalk(:,:,k).*XkF);
        FXkP = Y_abs_noisy(:,:,k).*exp(i*angle(FXk)); 
        XkF  = ifft2(FXkP)./Diagonalk(:,:,k);
    end
    %% PoPf
    XkOF = zeros(N11,N22);
    XkOF(N1start:N1end,N2start:N2end) = projection(XkF(N1start:N1end,N2start:N2end),imagetype);
    %% updated Xk
    Xk = Xk - XkO - XkF + 2*XkOF;
    %% Imagek
    Imagek  = zeros(N11,N22);
    Imagek(N1start:N1end,N2start:N2end) = projection(Xk(N1start:N1end,N2start:N2end),imagetype);
    %% update the mask with Xk
    for k = indexstart : P
        FXk         = fft2(Diagonalk(:,:,k).*Xk);
        FXkP        = Y_abs_noisy(:,:,k).*exp(i*angle(FXk)); 
        %% ONLY ANGLE
        DiagonalkP  = ifft2(FXkP) ./ exp(i*angle(Xk));
        %% fit the mask
        AngleInk  = angle(DiagonalkP(N1start:N1end,N2start:N2end));
        Diff      = abs(AngleInk - Angle0(N1start:N1end,N2start:N2end,k));
        Fit1      = (min(Diff,2*pi-Diff) <= anglerange);
        NoFit1    = 1 - Fit1;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if maskupdate == 1
            AngleUnfit = Angle0(N1start:N1end,N2start:N2end,k);
        elseif maskupdate == 3
            AngleUnfit = zeros(N1end-N1start+1,N2end-N2start+1);
            for index1 = 1 : N1end-N1start+1
                for index2 = 1 : N2end-N2start+1
                    if Fit1(index1,index2) == 0
                        disminus = abs(exp(i*AngleInk(index1,index2))-exp(i*(Angle0(index1+N1start-1,index2+N2start-1,k)-anglerange)));
                        displus  = abs(exp(i*AngleInk(index1,index2))-exp(i*(Angle0(index1+N1start-1,index2+N2start-1,k)+anglerange)));
                        if disminus <= displus
                            AngleUnfit(index1,index2) = Angle0(index1+N1start-1,index2+N2start-1,k)-anglerange;
                        else
                            AngleUnfit(index1,index2) = Angle0(index1+N1start-1,index2+N2start-1,k)+anglerange;
                        end
                    end
                end
            end
        end
        AngleInk   = AngleInk.*Fit1 + AngleUnfit.*NoFit1;
        Diagonalk(N1start:N1end,N2start:N2end,k) = exp(i*AngleInk); 
    end
    %% residual
    reres = 0;
    for k = 1 : P
        reres = reres + norm(abs(fft2(Diagonalk(:,:,k).*Imagek))-Y_abs_noisy(:,:,k),'fro')^2;
    end
    reres = sqrt(reres)/Y_abs_fro;
    iterres(iter)  = reres;
    %% display the results
    if ifdisplay == 1
        fprintf('%6.0f res = %6.4f  \n',iter,reres)
    end
    %% To be adaptive, terminate DRER if residual keeps increasing in five
    %% steps
    if ifadaptive == 1 
        if (imagetype==1) || (imagetype==2)
            %% step 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
% % %     if norm(Xk-Xpre,'fro')/norm(Xpre,'fro') < 0.01
% % %         break;
% % %     end
end

hiostep  = length(iterres);
increase = 0;


%% AER iteration
for iter = 1 : ifer
    %%
    Xpre = Xk;
    %% update X
    XkP   = Xk;
    for k = 1 : P
        FXk  = fft2(Diagonalk(:,:,k).*XkP);
        FXkP = Y_abs_noisy(:,:,k).*exp(i*angle(FXk)); %*(FXk./abs(FXk));
        XkP  = ifft2(FXkP)./Diagonalk(:,:,k);
    end
    %%%% update Xk to fit the the space domain constraint   
    Xk = zeros(N11,N22);
    Xk(N1start:N1end,N2start:N2end) = projection(XkP(N1start:N1end,N2start:N2end),imagetype);
    %% update the mask with Xk
    for k = indexstart : P
        FXk         = fft2(Diagonalk(:,:,k).*Xk);
        FXkP        = Y_abs_noisy(:,:,k).*exp(i*angle(FXk)); 
        %% ONLY ANGLE
        DiagonalkP  = ifft2(FXkP) ./ exp(i*angle(Xk));
        %% fit the mask
        AngleInk  = angle(DiagonalkP(N1start:N1end,N2start:N2end));
        Diff      = abs(AngleInk - Angle0(N1start:N1end,N2start:N2end,k));
        Fit1      = (min(Diff,2*pi-Diff) <= anglerange);
        NoFit1    = 1 - Fit1;
        if maskupdate == 1
            AngleUnfit = Angle0(N1start:N1end,N2start:N2end,k);
        elseif maskupdate == 3
            AngleUnfit = zeros(N1end-N1start+1,N2end-N2start+1);
            for index1 = 1 : N1end-N1start+1
                for index2 = 1 : N2end-N2start+1
                    if Fit1(index1,index2) == 0
                        disminus = abs(exp(i*AngleInk(index1,index2))-exp(i*(Angle0(index1+N1start-1,index2+N2start-1,k)-anglerange)));
                        displus  = abs(exp(i*AngleInk(index1,index2))-exp(i*(Angle0(index1+N1start-1,index2+N2start-1,k)+anglerange)));
                        if disminus <= displus
                            AngleUnfit(index1,index2) = Angle0(index1+N1start-1,index2+N2start-1,k)-anglerange;
                        else
                            AngleUnfit(index1,index2) = Angle0(index1+N1start-1,index2+N2start-1,k)+anglerange;
                        end
                    end
                end
            end
        end
        AngleInk   = AngleInk.*Fit1 + AngleUnfit.*NoFit1;
        Diagonalk(N1start:N1end,N2start:N2end,k) = exp(i*AngleInk); 
    end 
    %% residual
    reres = 0;
    for k = 1 : P
        reres = reres + norm(abs(fft2(Diagonalk(:,:,k).*Xk))-Y_abs_noisy(:,:,k),'fro')^2;
    end
    reres = sqrt(reres)/Y_abs_fro;
    iterres(hiostep+iter)  = reres;
    %% display the results
    if ifdisplay == 1
        fprintf('%6.0f res = %6.4f \n',iter,reres)
    end
    if ifadaptive == 1
        if (norm(Xk-Xpre,'fro')/norm(Xpre,'fro') < 0.0005)
            break;
        end
    end
end

erstep  = length(iterres) - hiostep;
Xrec    = Xk;





