% 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.5 for iterative 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, L=2l+1)  
               

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

   % number of iterations in algorithm 2.5 depends on the reduction
   % factor 
   if R==2 iter=12; 
   elseif R==3 iter=40; 
   elseif R==4 iter=100; 
   elseif R==5 iter=50; 
   elseif R==6 iter=90; 
   end
                                                                       
[N,M,Nc]=size(Y)
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')


%illustration of optimal coil sensitivities from the data 
% Figure 5 in the paper
mrec1=mrec0;
for k1=1:N 
    for k2=1:M 
        if mrec0(k1,k2)>0.005 mrec1(k1,k2)=1/mrec0(k1,k2);
        else mrec1(k1,k2)=0;
        end
    end
end
%optimal coil sensitivities
sopt=zeros(N,M,Nc); sopthat=zeros(N,M,Nc);sopthats=zeros(30,30,Nc);
for j=1:Nc
sopt(:,:,j)=Ycheck(:,:,j).*mrec1;
sopthat(:,:,j)= fftshift(fft2(fftshift(sopt(:,:,j))))./sqrt(M*N);
sopthats(:,:,j)=sopthat(floor(N/2)-15:floor(N/2)+14,floor(M/2)-15:floor(M/2)+14,j);
end


figure('name','optimal sensitivities'); 
imshow3(abs(sopt),[],[1,Nc]); axis off;
set(gca,'FontSize',10);
%colormap(gray(256)); 
colormap default
colorbar
%saveas(gcf,'optimal_sensitivities.eps')

sopthats=sopthats./max(abs(sopthats(:)));
figure('name','optimal sensitivities in k-space'); 
imshow3(abs(sopthats),[],[1,Nc]); axis off;
set(gca,'FontSize',10);
colormap default
colorbar
%saveas(gcf,'optimal_sensi_k.eps')



%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.3 in the paper 
% using the iterative Algorithm 2.5

[mrec,Yrec] = multichannel_iterative_mri1(YP,P,M1,L,beta,iter,R);
 mrec=abs(mrec);
 mrec=mrec./max(mrec(:)); 
 
 % Figure of reconstructed image mrec (1% of largest values projected to max of the other
 % values)
 mrec1= scaling_for_illustration(mrec,100);
 figure('name','Reconstructed image using the iterative algorithm'),
 imagesc(mrec1); axis off;
 set(gca,'FontSize',20);
 colormap(gray);
 set(gca,'LooseInset',get(gca,'TightInset'));
 %saveas(gcf,'MOCCA_R2_L3_it12.eps')
 
% Comparison with reconstruction from complete data
 m1=mrec0./max(mrec0(:));    % normalization of the reference image (full data) on [0,1]
 p = psnr(m1,abs(mrec))           % PSNR value of the error
 ssimval = ssim(m1,abs(mrec))     % SSIM of the error
 
 %Error map
 error=abs(mrec-m1);  % relative error (since max(abs(m1))=1)

 figure('name','error map'), imagesc(error); axis off;
 colormap default
 colorbar %('off')
 set(gca,'FontSize',18);
 %saveas(gcf,'error_R2_M8.eps')
 %exportgraphics(gca,'error_R2_L3_it12bar.eps','Resolution',600)
 
