function [mrec,Yrec] = multichannel_iterative_mri1(YP,P,M1,L,beta,iter,R)
% function  multichannel algorithm (Algorithm 2.3 and Algorithm 2.5 in 
%  multichannel  paper) and SoS condition
%  for reconstruction of MRI images from incomplete k-space data
%
%
%  Input: YP acquired coil image data YP=(Y^1.*P ...Y^Nc.*P ) of size
%          NxMxNc, NxM dimension of the image to reconstruct (N, M even)
%         P projection NxM matrix, contains ones if corresponding Y^j entry is
%         acquired, zeros otherwise
%         Nc = number of coils
%         L size of support of coil sensitivities in k-space, L odd
%         M1 = calibration matrix dimension in calibration area providing a
%         M1+L-1 box (M1 even)
%         N x M dimension of the image to reconstruct, N even! M even!
%         iter upper bound for number of iterations in Algorithm 2.5 to solve the least
%         squares problem
%         beta regularization parameter in the least squares problem 
%         R reduction factor (for unacquired data)
%
% Output: mrec reconstructed image of size NxM,
%         Yrec =(Yrec^1, ... , Yrec^Nc) of size NxM*Nc reconstructed k-space
%         data


[N,M,Nc]=size(YP);

%Build Mocca matrix
AMocca=Moccamatrix(YP,M1,L);

% compute the normalized sensitivities using the calibration data
srec=multichannel_coils_recon(AMocca,L,N,M,Nc,R); 

%solve the least squares system in Algorithm 2.3 iteratively
%with Algorithm 2.5

%iteration initialization
m0=zeros(N,M); m1=m0;  Yrec=zeros(N,M,Nc);

% initial approximation of m from incomplete k-space data YP

for j=1:Nc
    cY=sqrt(N*M).*ifftshift(ifft2(ifftshift(YP(:,:,j))));
    m1=m1+ conj(srec(:,:,j)).*cY;
end

%iteration step (Yrec and m are updated at each step)
ind=0;
while and(abs(m0(:)-m1(:))>10^(-10),ind<iter)
    ind=ind+1;
    m2=zeros(N,M);
    for j=1:Nc
        cY=(fftshift(fft2(fftshift(m1.*srec(:,:,j)))))./sqrt(M*N);
        Yrec(:,:,j)=(1-beta).*(ones(N,M)-P).*cY+ YP(:,:,j);
        m2= m2+ sqrt(M*N).*conj(srec(:,:,j)).*ifftshift(ifft2(ifftshift(Yrec(:,:,j))));
    end
    m0=m1; m1=m2;
end
ind=ind

% illustrate coil sensitivities

mrec=m1;
index=find(abs(mrec) >=0.004);
signmrec=zeros(N*M,1);
mrec1=reshape(mrec,N*M,1);
signmrec(index)=sign(mrec1(index));
signmrec=reshape(signmrec,N,M);
for j=1:Nc srec(:,:,j)=srec(:,:,j).*signmrec; end

%illustrate coils sensitivities

figure('name','magnitude of coil sensitivities after multiplication with sign(m)') 
imshow3(abs(srec),[],[2,Nc/2]) 
axis off
set(gca,'FontSize',20)
colormap default
colorbar
set(gca,'FontSize',20);
%saveas(gca,'coilabs1_L=5_m.eps','epsc')
%colormap default

figure('name','coil sensitivities (phase angle) after multiplication with sign(m)'), 
imshow3(angle(srec),[],[2,Nc/2]) 
axis off
set(gca,'FontSize',20)
colorbar
phasemap
set(gca,'FontSize',20);
%saveas(gca,'coilphase1_L=5_m.eps','epsc')
end

