function [mrec,Yrec] = multichannel_fast_mri1(YP,P,P1,M1,L,Nc,beta,R)
% reconstruction from incomplete MRI data (Algorithm 2.3)
%    using the direct solution method in Algorithm 2.7
%
%  Input: YP  acquired coil image data YP=(Y^1.*P ...Y^Nc.*P ) of size N times M*Nc
%         P   projection matrix including the calibration area
%         P1  projection matrix, contains ones according to the structure fixed by
%         L= 2l+1 (odd), length of support of sensitivities in k-space
%         M1 = calibration matrix dimension in calibration area providing a
%         M1+L-1 box (M1 even)
%         Nc = number of coils
%         R reduction factor 
%         depending on R (resp. the projection matrix)we need M,N even
%         for R=2 N,M being multiples of 2, 
%         for R=3 we need N being a multiple of 2 and M being a multiple of 6
%         for R=5 (each second row and each second column acquired) N,M being multiples of 4,
%         for R=4 (each fourth colun acquired) we need N to be a multiple of 2 and  M being
%         multiple of 8, 
%         for R=6 we need N being a multiple of 4 and M being a multiple of 6
% Output: reconstructed image mrec, reconstructed data Yrec

[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); 



%Compute the equation system in step 5 of Algorithm 2.3 directly using the
%special structure of the projection matrix P1 

%righthand side
Rs=zeros(N,M);
for j=1:Nc
Rs=Rs+conj(srec(:,:,j)).*(ifftshift(ifft2(ifftshift(P1.*YP(:,:,j)))))*sqrt(N*M);
end
% Define G matrix (which decomposes into smaller matrices according to the structure of
% P1 and solve the equation system to compute m
                   
if R==2  %each second column acquired
% Define G matrix (which decomposes into N*M/2 (2x2)-matrices) and solve the system to compute m
m0=zeros(N,M);
for k1=1:N
    for k2=1:M/2
        % define 2x2 partial matrices of G
        G=zeros(2,2); re=zeros(2,1); %initialize G matrix and right-hand side
        for j=1:Nc
            s=[srec(k1,k2,j) srec(k1,k2+M/2,j)];
            G=G+s'*s;
        end
        G= G+ 2*beta*eye(2)/(N*M);
        re=2*[Rs(k1,k2); Rs(k1,k2+M/2)]; %right hand side
        %solve 2x2 system
        mm=linsolve(G,re);
        m0(k1,k2)=mm(1);
        m0(k1,k2+M/2)=mm(2);
    end
end
elseif R==3  %each third column acquired
% Define G matrix (which decomposes into N*M/3 (3x3)-matrices) and solve the system to compute m
m0=zeros(N,M);
for k1=1:N
    for k2=1:M/3
        % define 2x2 partial matrices of G
        G=zeros(3,3); re=zeros(3,1); %initialize G matrix and right-hand side
        for j=1:Nc
            s=[srec(k1,k2,j) srec(k1,k2+M/3,j) srec(k1,k2+2*M/3,j)];
            G=G+s'*s;
        end
        G= G+ 3*beta*eye(3)/(N*M);
        re=3*[Rs(k1,k2); Rs(k1,k2+M/3); Rs(k1,k2+2*M/3)]; %right hand side
        %solve 2x2 system
        mm=linsolve(G,re);
        m0(k1,k2)=mm(1);
        m0(k1,k2+M/3)=mm(2);
        m0(k1,k2+2*M/3)=mm(3);
    end
end
elseif R==4        % each fourth column acquired
                   % define G matrix (which decomposes into N*M/4 (4x4)-matrices) 
                   % and solve the system to compute m
m0=zeros(N,M);
for k1=1:N
    for k2=1:M/4
        % define 4x4 partial matrices of G
        G=zeros(4,4); re=zeros(4,1); %initialize G matrix and right-hand side
        for j=1:Nc
            s=[srec(k1,k2,j) srec(k1,k2+M/4,j) srec(k1,k2+M/2,j) srec(k1,k2+3*M/4,j)];
            G=G+s'*s;
        end
        G= G+ 4*beta*eye(4)/(N*M);
        re=4*[Rs(k1,k2); Rs(k1,k2+M/4); Rs(k1,k2+M/2);  Rs(k1,k2+3*M/4)]; %right hand side
        %solve 4x4 system
        mm=linsolve(G,re);
        m0(k1,k2)=mm(1);
        m0(k1,k2+M/4)=mm(2);
        m0(k1,k2+M/2)=mm(3);
        m0(k1,k2+3*M/4)=mm(4);
    end
end
elseif R==5        % each second row and each second column acquired
                   % define G matrix (which decomposes into N*M/4 (4x4)-matrices) 
                   % and solve the system to compute m
m0=zeros(N,M);
for k1=1:N/2
    for k2=1:M/2
        % define 4x4 partial matrices of G
        G=zeros(4,4); re=zeros(4,1); %initialize G matrix and right-hand side
        for j=1:Nc
            s=[srec(k1+N/2,k2,j) srec(k1+N/2,k2+M/2,j) srec(k1,k2,j)  srec(k1,k2+M/2,j)];
            G=G+s'*s;
        end
        G= G+ 4*beta*eye(4)/(N*M);
        re=4*[Rs(k1+N/2,k2); Rs(k1+N/2,k2+M/2); Rs(k1,k2);  Rs(k1,k2+M/2)]; %right hand side
        % solve 4x4 system
        %mm=linsolve(G,re);
        mm=pinv(G)*re;
        m0(k1+N/2,k2)=mm(1);
        m0(k1+N/2,k2+M/2)=mm(2);
        m0(k1,k2)=mm(3);
        m0(k1,k2+M/2)=mm(4);
    end
end
elseif R==6  % each second row and each third column acquired
             % define G matrix (which decomposes into N*M/6 (6x6)-matrices) and solve the system to compute m
m0=zeros(N,M); %condition=zeros(N/2,M/2);
for k1=1:N/2
    for k2=1:M/3
        %construct G matrix of size 6 x 6
        G=zeros(6,6); re=zeros(6,1); %initialize G matrix and right-hand side
        for j=1:Nc
            s=[srec(k1+N/2,k2,j) srec(k1+N/2,k2+M/3,j) srec(k1+N/2,k2+2*M/3,j) srec(k1,k2,j)  srec(k1,k2+M/3,j) srec(k1,k2+2*M/3,j)];
            G=G+s'*s;
        end
        G= G+ 6*beta*eye(6)/(N*M);
        re=6*[Rs(k1+N/2,k2); Rs(k1+N/2,k2+M/3);Rs(k1+N/2,k2+2*M/3);  Rs(k1,k2);  Rs(k1,k2+M/3); Rs(k1,k2+2*M/3)]; %right hand side
        %solve 6x6 system
        mm=linsolve(G,re);
        %mm=pinv(G)*re;
        m0(k1+N/2,k2)=mm(1);
        m0(k1+N/2,k2+M/3)=mm(2);
        m0(k1+N/2,k2+2*M/3)=mm(3);
        m0(k1,k2)=mm(4);
        m0(k1,k2+M/3)=mm(5);
        m0(k1,k2+2*M/3)=mm(6);
    end
end
end

%renormalization of m not needed
%mrec=m0;


% Improvement by one iteration step  of Algorithmus 2.5
Yrec=zeros(N,M,Nc);
for k=1:1
    m1=zeros(N,M);
for j=1:Nc
        cY=fftshift(fft2(fftshift(m0.*srec(:,:,j))))./sqrt(M*N); 
        Yrec(:,:,j)=(1-beta).*(ones(N,M)-P).*cY+ YP(:,:,j);
        m1= m1+ sqrt(M*N).*conj(srec(:,:,j)).*ifftshift(ifft2(ifftshift(Yrec(:,:,j))));
end
m0=m1;
end
mrec=m0;
end

