%% Script to reproduce Figure 4

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                        
% Nearfield ACoustic HOlography with Sparse regularization (NACHOS)
% Version 1.0                                    
%     
% Copyright 2012 Antoine Peillot, Fran蔞is Ollivier, Gilles Chardon, 
%                Laurent Daudet, Nancy Bertin, R幦i Gribonval
% 
% For all details please refer to README.TXT
%
% This software is a free software distributed under the terms of the GNU 
% Public License version 3 (http://www.gnu.org/licenses/gpl.txt). You can 
% redistribute it and/or modify it under the terms of this licence.                                        
%                                                                                                                                                      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% load data

Fmax = 4000;

SPECPVibroD.SelectedFreqIndices = [327 981 2588 5582];
SPECPVibroD.RelevantFreqIndices = [];

Nbx = size(VibroDANTP.Mask,1);
Nby = size(VibroDANTP.Mask,2);

% Algorithm & target reconstruction specifications
ALGOP.Name = 'l1_constr_spgl1';
ALGOP.StoppingCriterion = 27;
ALGOP.DicoSizeX = 2*Nbx;
ALGOP.DicoSizeY = 2*Nby;
ALGOP.StepX = 0.0085;
ALGOP.StepY = 0.0082;
ALGOP.DicoXLength = ALGOP.DicoSizeX*ALGOP.StepX;
ALGOP.DicoYLength = ALGOP.DicoSizeY*ALGOP.StepY;

ALGOP.Mask = VibroDANTP.Mask;
DisplayMask = ALGOP.Mask;
DisplayMask(DisplayMask==0) = NaN;

% prepare loops
NbMethods = 6;

ODSToPlot = cell(NbMethods,length(SPECPVibroD.SelectedFreqIndices));

Correlations = zeros(NbMethods,length(SPECPVibroD.SelectedFreqIndices));

%% Loop on modes to be imaged
for ii = 1:length(SPECPVibroD.SelectedFreqIndices)
    
    % Prepare groundtruth image

    ODSFreqIdx = SPECPVibroD.SelectedFreqIndices(ii);
    ODS        = reshape(SpectraVibroD(:,ODSFreqIdx),Nby,Nbx).';
    
    Correlations(1,ii) = 1;
    
    % Reconstruction with Tikhonov reg. from regular antenna measurements
    nfft = 128; % the fft2d is 128x128 
    DownSampleParam = [1,1,48,1,1,40];
    [Spectrum, REGANTP1920] = select_datasubset(SpectraRegAntD, RegDANTP, ODSFreqIdx,0,'none',DownSampleParam);
    ODS_RegAnt1920_Tikhonov  = apply_tikhonov_reconstruction(Spectrum,SPECPRegAntD,ODSFreqIdx,REGANTP1920,VibroDANTP,nfft);
    Correlations(2,ii) = compute_correlation_2d(ODS,ODS_RegAnt1920_Tikhonov.*ALGOP.Mask);
    
    DownSampleParam = [3,4,47,3,4,39];
    [Spectrum, REGANTP120] = select_datasubset(SpectraRegAntD, RegDANTP, ODSFreqIdx,0,'regular',DownSampleParam);
    ODS_RegAnt120_Tikhonov  = apply_tikhonov_reconstruction(Spectrum,SPECPRegAntD,ODSFreqIdx,REGANTP120,VibroDANTP,nfft);
    Correlations(3,ii) = compute_correlation_2d(ODS,ODS_RegAnt120_Tikhonov);
    
    % Sparse Reconstructions from regular antenna measurements (120 micros over and out of the plate)
    ODS_RegAnt120_BPDN = apply_sparse_reconstruction(Spectrum,SPECPRegAntD,ODSFreqIdx,REGANTP120, ALGOP,VibroDANTP);
    Correlations(4,ii) = compute_correlation_2d(ODS,ODS_RegAnt120_BPDN.*ALGOP.Mask);
    
    % Sparse Reconstructions from regular antenna measurements (120 micros
    % over the plate)
    DownSampleParam = [7,3,40,3,3,32];
    [Spectrum, REGANTP120p] = select_datasubset(SpectraRegAntD, RegDANTP, ODSFreqIdx,0,'regular',DownSampleParam);
    ODS_RegAnt120p_BPDN = apply_sparse_reconstruction(Spectrum,SPECPRegAntD,ODSFreqIdx,REGANTP120p, ALGOP,VibroDANTP);
    Correlations(5,ii) = compute_correlation_2d(ODS,ODS_RegAnt120p_BPDN.*ALGOP.Mask);
    
    % Reconstruction with BPDN, from rand antenna measurements
    DownSampleParam = 120;
    [Spectrum,RANDANTP120] = select_datasubset(SpectraRegAntD, RegDANTP,ODSFreqIdx,1,'random',DownSampleParam); % TODO 
    ODS_RandAnt120_BPDN = apply_sparse_reconstruction(Spectrum,SPECPRegAntD,ODSFreqIdx,RANDANTP120, ALGOP,VibroDANTP);
    Correlations(6,ii) = compute_correlation_2d(ODS,ODS_RandAnt120_BPDN.*ALGOP.Mask);
           
    % store for display
    ODSToPlot{1,ii} = rot90(Cpx2Surf(ODS));
    ODSToPlot{2,ii} = rot90(ODS_RegAnt1920_Tikhonov);
    ODSToPlot{3,ii} = rot90(ODS_RegAnt120_Tikhonov);
    ODSToPlot{4,ii} = rot90(ODS_RegAnt120_BPDN);
    ODSToPlot{5,ii} = rot90(ODS_RegAnt120p_BPDN);
    ODSToPlot{6,ii} = rot90(ODS_RandAnt120_BPDN);   
    
end

% create figure 4
f4=figure(4);

set(f4,'OuterPosition',[500, 500, 1000, 800]);
plot_ODS(ODSToPlot,Correlations,rot90(DisplayMask));

subplot(length(SPECPVibro.SelectedFreqIndices),NbMethods,1);
title('(a) Groundtruth');
subplot(length(SPECPVibro.SelectedFreqIndices),NbMethods,2);
title('(b) TikhREG1920');
subplot(length(SPECPVibro.SelectedFreqIndices),NbMethods,3);
title('(c) TikhREG120');
subplot(length(SPECPVibro.SelectedFreqIndices),NbMethods,4);
title('(d) L1REG120');
subplot(length(SPECPVibro.SelectedFreqIndices),NbMethods,5);
title('(e) L1REG120in');
subplot(length(SPECPVibro.SelectedFreqIndices),NbMethods,6);
title('(f) L1RANfREG120in');

subplot(length(SPECPVibro.SelectedFreqIndices),NbMethods,1);
ylabel([num2str(round(SPECPVibro.FreqValues(SPECPVibro.SelectedFreqIndices(1)))),' Hz']);
subplot(length(SPECPVibro.SelectedFreqIndices),NbMethods,7);
ylabel([num2str(round(SPECPVibro.FreqValues(SPECPVibro.SelectedFreqIndices(2)))),' Hz']);
subplot(length(SPECPVibro.SelectedFreqIndices),NbMethods,13);
ylabel([num2str(round(SPECPVibro.FreqValues(SPECPVibro.SelectedFreqIndices(3)))),' Hz']);
subplot(length(SPECPVibro.SelectedFreqIndices),NbMethods,19);
ylabel([num2str(round(SPECPVibro.FreqValues(SPECPVibro.SelectedFreqIndices(4)))),' Hz']);