% test of Algorithm 2.3 in the paper 'Model-Based Multi-Channel Blind Deconvolution from Incomplete Fourier 
% Data' (by G. Plonka and Y. Riebe)
% here employing Algorithm 2.7 for fast solution of
% the equation system in step 5 of Algorithm 2.3. 


% Output: -8 coil images (inverse Fourier transform from complete coil k-space data)
%         -reconstructed image mrec0 from complete k-space data computed with SOS condition
%         -reconstruced image mrec from incomplete k-space data
%         -pointwise error compared to mrec0
clc
clear all
close all


%Choice of parameters

M1=20;          % fixed size of calibration area: (M1+l-1)x(M1+L-1) box around center
beta=0.0001;    % fixed regularization parameter in iteration of Algorithm 2.3
R=2;            % Reduction factor for incomplete data (1=complete data, 2=each second column, 
                % 3= each third column, 4= each fourth column, 5=each second row and each second column, 
                % 6=each second row and each third column acquired)
L=5;            % support size of coil sensitivities in k-space (odd, larger than 5)  
               

                
 load brain_8ch.mat;   % data set (complete)
 DATA = DATA/max(max(max(abs(ifftshift(ifft2(ifftshift(DATA))))))) + eps;
 Y = DATA./max(abs(DATA(:)));

 % it is assumed that N and M are suitable lengths
[N,M,Nc]=size(Y)
if R==3 Y=Y(:,1:3*floor(M/3),:); [N,M,Nc]=size(Y); end
if R==3 Y=Y(:,1:3*floor(M/3),:); [N,M,Nc]=size(Y); end
Ycheck=sqrt(M*N)*ifftshift(ifft2(ifftshift(Y)));
    
% illustration of coil images (with proper scaling)
% see Figure 4 in the paper
Ycheck1= scaling_for_illustration(Ycheck,150);
figure('name','coil images'); 
imshow3(abs(Ycheck1),[],[1,Nc]); axis off;
set(gca,'FontSize',20);
colormap(gray(256)); 
set(gca,'FontSize',20);

%reconstruction mrec0 from complete k-space data with sos condition
mrec0=zeros(N,M);         % initialize mrec0
for j=1:Nc
    Yz=Ycheck(:,:,j);
    mrec0=mrec0+conj(Yz).*Yz;   
end
mrec0=sqrt(mrec0); 

% figure of reconstructed image mrec0 from complete k-space data (with SoS)
mrec00= scaling_for_illustration(mrec0,100);
figure('name','image reconstruction from full data with SOS'),
imagesc(mrec00); axis off;
set(gca,'FontSize',10);
colormap(gray);
%saveas(gcf,'SOS_reconstruction.png')


%input projection matrix (zero for each unacquired entry)
[P,P1]= projection(N,M,R,M1,L);

% Determine given (incomplete) k-space data
YP = Y;
for j=1:Nc
YP(:,:,j) = Y(:,:,j).*P;  % YP only contains the acquired k-space data
end      


% Multichannel reconstruction according to Algorithm 2.5 in the paper using
% fast Algorithm 2.7

[mrec,Yrec] = multichannel_fast_mri1(YP,P,P1,M1,L,Nc,beta,R);
mrec=abs(mrec);


% Comparison with reconstruction from complete data
 m1=mrec0./max(mrec0(:));    % normalization of the reference image (full data) on [0,1]
 mrec=mrec./max(mrec(:));    % normalization on [0,1]
 p = psnr(m1,mrec)           % PSNR value of the error
 ssimval = ssim(m1,mrec)     % SSIM of the error


 
 %Error map
 error=abs(mrec-m1)/max(m1(:));  %relative error
 figure('name','error map'), imagesc(error); axis off;
 colormap default
 colorbar
 set(gca,'FontSize',20);
 %caxis([0 0.2]);
 
 % Figure of reconstructed image mrec (1% of largest values projected to max of the other
 % values)
 mrec=abs(mrec);
 mrec=mrec./max(mrec(:));
 mrec1= scaling_for_illustration(mrec,100);
 
 
 figure('name','reconstructed image with Mocca using direct computation'),
 imagesc(mrec1); axis off;
 set(gca,'FontSize',20);
 colormap(gray);
 set(gca,'LooseInset',get(gca,'TightInset'));
 %saveas(gcf,'MOCCA2_R4_L7_Ns47.eps')
 

 
