function BlindDecSelectAIF(working_dir,form_name,offline)

working_dir
form_name



addpath('./PerfLabUtils');
addpath('./BlindDec');
addpath('./NonBlindDec');
addpath('./General');


NextNextstep=NextXmlFile(form_name,2);
PBfile = fullfile(working_dir,'progressBlindDec.txt');

%% Read GUI parameters
fileName=fullfile(working_dir,[form_name '.dat']); %0175_SequenceViewer.dat'
objekty =xml2PerfLab(fileName);

aif_estim = FromListObj(objekty,'aif_estim'); %% AIF estimation method
BlindModel = FromListObj(objekty,'BlindModel'); %% IRF model for blind dec.
AIFModel = FromListObj(objekty,'AIFModel'); %% AIF model for blind dec.

NumBBB = str2num(objekty.NumBBB.value.Text);

%% AIF filename - for loading
%AIFindex= str2num( objekty.AIFsList.options.option.Text) % index of the selected AIF in the list
AIFindex = str2num( FromListObj(objekty,'AIFsList') );


%% input ROI string
input_string_ROI=objekty.sv.value.Text	% string defining the selected ROIs

%% read polygon / voxel coordinates from the returned string
%zvetseni = 4;
Polygons = ReadROIVoxelString(input_string_ROI);

%% Input data file
InpFile = fullfile(working_dir,'conv','data.mat')

%% Output dir for mat file
imgfolder = 'aif'
outpath=fullfile(working_dir,imgfolder)
if ~exist(outpath,'dir');
    mkdir(outpath)
end
OutFile = fullfile(outpath,'data.mat')

%% Output plot name
OutputPlotName = fullfile(working_dir,'graph_aif.svg');

%% Output dir with parametric maps for viewer
imgfolder = 'viewer_par'
outpath=fullfile(working_dir,imgfolder)
delete([outpath '/*.*']);
if ~exist(outpath,'dir');
    mkdir(outpath)
end
ViewerParDir = outpath

%% output dir for parametric maps in Matlab format
MapFile = fullfile(outpath,'maps.mat')


%% Size of SequenceViewer image (for scaling)
fileName=fullfile(fileparts(pwd),[form_name '.xml'])
objekty =xml2PerfLab(fileName); %OMmod
svX=objekty.sv.width.Text
svY=objekty.sv.height.Text
svX=str2num(svX)
svY=str2num(svY)

load(InpFile);


switch aif_estim
    
    case 'load'
        %% Load AIF from a file
        
        %% Need for temporal upsampling?
        if info.acq.Ts > 1.5
            disp('Upsampling to Ts=1s')
            t = info.acq.Ts *(0:size(data2,2)-1) / 60; %min
            t_new = 0:1/60:max(t); %min (temp. sampling 1 sec)
            
            d2 = DataCells2Mat(data2);
            Zs = size(d2,3);
            Ys = size(d2,1);
            Xs = size(d2,2);
            d3 = zeros(Ys,Xs,Zs,length(t_new));
            for z=1:Zs
                for y=1:Ys
                    for x=1:Xs
                        d3(y,x,z,:) = interp1(squeeze(t),squeeze(d2(y,x,z,:)),t_new,'spline');
                    end
                end
            end
            data2 = Mat2DataCells(d3);
            info.acq.Ts = 1;
            movefile(InpFile, [InpFile(1:end-4) '_oldTs.mat']);
            save(InpFile,'data1','data2','info','-v7.3');
        end
        
        
        %% temporary file for read2matlab
        tmpfile = fullfile(working_dir,'iitmp.asci')
        
        %AIFsList = AIFsList(3:end)
        %AIFsList = strrep(AIFsList,'#','/');
        %AIFsList = strrep(AIFsList,'@','.');
        aiflist = fopen([working_dir '/aiflist.txt'],'r');
        for line=1:AIFindex
            AIFsList = fgetl(aiflist)
        end
        fclose(aiflist);
        
        
        [pathstr, name, ext] = fileparts(AIFsList)
        
        if strcmp(ext,'.real')
            % XITE format
            ok = 1;
            aif0 = [];
            frame = 1;
            if exist('offline')
                global AifsDir
                aifnmg = AIFsList;
            else
                aifnmg = [working_dir '../../../../_storage/AIFs/' AIFsList];
            end
            while ok
                try
                    aifnm = [aifnmg ':' int2str(frame)];
                    aif0 = [aif0 real2matlab(aifnm,tmpfile)];
                catch
                    ok = 0;
                end
                frame = frame + 1;
                
                %                 if frame>length(data2)
                %                     break
                %                 end
            end
            aif = cell(1,3);
            aif{1} = aif0 / sum(aif0); % has to be normalized for correct scaling below
        else
            % Matlab format
            if exist('offline')
                global AifsDir
                %aiffile = load(fullfile(AifsDir, AIFsList));
                aiffile = load(AIFsList);
            else
                aiffile = load([working_dir '../../../../_storage/AIFs/' AIFsList]);
            end
            aif = aiffile.aif;
            aif{1} = aif{1}(:)';
            
            [Cm, indm] = max(aif{1});
            sprintf('Maximum of the imported AIF at index location: %d', indm)
            sprintf('Maximum of the imported AIF at time location : %f s', indm*aiffile.info.acq.Ts)
            
            %% resample and impose the same length
            t =             info.acq.Ts *(0:size(data2,2)-1) / 60;
            t_aif = aiffile.info.acq.Ts *(0:length(aif{1})-1) / 60;
            
            if (max(t)~=max(t_aif)) || (length(t)~=length(t_aif))
                % interpolate
                t=t(t<max(t_aif)); %!!!
                max(t)
                max(t_aif)
                %save(OutFile,'t_aif','aif','t');
                try
                    aif_res = interp1(t_aif,aif{1},t,'spline');
                    disp('interpolation done')
                catch
                    aif_res = interp1(t_aif,aif{1},t(t<max(t_aif)));
                    deltaL = length(t)-length(aif_res);
                    aif_res = [aif_res(:); zeros(deltaL,1)]';
                end
                aif{1} = aif_res;
            end
        end
        
        % data2 = cell(1,1);
        %data2{1} = data3{100};
        
        % save
        %save([DataFile1 filename_aif], 'data1','data2','info', 'aif');
        
        info.aif.method = 'imported';
        info.aif.model = 'unknown';
        %info.aif.scaling = ScalingMethod;
        
        tissue = aif;
        
        %% Plot AIF for next step:
        t = info.acq.Ts *(0:size(data2,2)-1) / 60;
        h=figure;
        size(t)
        size(aif{1})
        if length(aif{1}) < length(t)
            aif{1} = [aif{1} zeros(1,length(t)-length(aif{1}))];
        end
        if length(aif{1}) > length(t)
            aif{1} = aif{1}(1:length(t));
        end
        plot(t,aif{1}); xlabel('time [min]'); ylabel('AIF [a.u.]')
        set(h,'position',[100 100 900 400])
        %print(h,'-dpng',OutputPlotName);
        print(h,'-dsvg',OutputPlotName);
        
        
    case 'meas'
        %% Measured AIF
        
        
        %% Need for temporal upsampling?
        %         if info.acq.Ts > 1.5
        %             disp('Upsampling to Ts=1s')
        %             t = info.acq.Ts *(0:size(data2,2)-1) / 60; %min
        %             t_new = 0:1/60:max(t); %min (temp. sampling 1 sec)
        %
        %             d2 = DataCells2Mat(data2);
        %             Zs = size(d2,3);
        %             Ys = size(d2,1);
        %             Xs = size(d2,2);
        %             d3 = zeros(Ys,Xs,Zs,length(t_new));
        %             for z=1:Zs
        %                 for y=1:Ys
        %                     for x=1:Xs
        %                         d3(y,x,z,:) = interp1(squeeze(t),squeeze(d2(y,x,z,:)),t_new,'spline');
        %                     end
        %                 end
        %             end
        %             data2 = Mat2DataCells(d3);
        %             info.acq.Ts = 1;
        %             movefile(InpFile, [InpFile(1:end-4) '_oldTs.mat']);
        %             save(InpFile,'data1','data2','info','-v7.3');
        %         end
        
        
        
        X = size(data2{1},2)
        Y = size(data2{1},1)
        Z = size(data2,1)
        
        %% Create masks and extract tissue curves to 'tissue' - for voxels
        tissue = {};
        frames = size(data2,2);
        NumPoly = size(Polygons,1);
        Polygons
        for ind=1:NumPoly
            %mask = poly2mask(Polygons{ind,1}/zvetseni,Polygons{ind,2}/zvetseni,size(data2{1},1),size(data2{1},2));
            ind
            Polygons{ind,1}
            Polygons{ind,2}
            %% FOR ROIs
            mask = zeros(Y,X,Z);
            mask(:,:,Polygons{ind,5}) = poly2mask(Polygons{ind,1}/svX*X,Polygons{ind,2}/svY*Y, Y,X);
            %% FOR POINTs
            %             mask = zeros(size(data2{1}));
            %             mask(round(Polygons{ind,2}(1)/svY*Y), round(Polygons{ind,1}(1)/svX*X))  = 1;
            tissue{ind,2} = logical(mask);
            tissue{ind,3} = [ Polygons{ind,1}/svX*X, Polygons{ind,2}/svY*Y ];
            curve = zeros(1,frames); %zeros(frames,1);
            for f=1:frames
                curve(f) = sum(sum(sum(mask .* DataCells2Mat(data2(:,f)) ))) / sum(sum(sum(mask)));
            end
            tissue{ind,1} = curve;
        end
        
        %% Take average of extracted curves as the AIF
        aif = zeros(size(tissue{1,1}));
        for ind = 1:size(tissue,1)
            aif = aif + tissue{ind,1};
        end
        aif = aif / size(tissue,1);
        aif = {aif};
        
        info.aif.method = 'select';
        info.aif.model = 'none';
        %info.aif.scaling = ScalingMethod;
        
        %% Plot AIF for next step:
        t = info.acq.Ts *(0:size(data2,2)-1) / 60;
        h=figure;
        size(t)
        size(aif{1})
        if length(aif{1}) < length(t)
            aif{1} = [aif{1} zeros(1,length(t)-length(aif{1}))];
        end
        if length(aif{1}) > length(t)
            aif{1} = aif{1}(1:length(t));
        end
        plot(t,aif{1}); xlabel('time [min]'); ylabel('AIF [a.u.]')
        set(h,'position',[100 100 900 400])
        %print(h,'-dpng',OutputPlotName);
        print(h,'-dsvg',OutputPlotName);
        
    case 'automeas'
        %% Measured AIF
        
        X = size(data2{1},2)
        Y = size(data2{1},1)
        Z = size(data2,1)
        
        %% Create masks and extract tissue curves to 'tissue' - for voxels
        tissue = {};
        frames = size(data2,2);
        NumPoly = size(Polygons,1);
        % %         Polygons
        mask = zeros(Y,X,Z); %mask = zeros(Y,X,Z);
        for ind=1:NumPoly
            %% FOR ROIs
            masktmp = zeros(Y,X,Z);
            masktmp(:,:,Polygons{ind,5}) = poly2mask(Polygons{ind,1}/svX*X,Polygons{ind,2}/svY*Y, Y,X);
            mask(masktmp>0) = 1;
        end
        %         data = DataCells2Mat(data2);
        %         data = permute(data,[4 1 2 3]);
        %         fullmask = permute(mask,[4 1 2 3]);
        %         data = data(logical(repmat(fullmask,[frames 1 1 1])));
        %         voxels = length(data)/frames;
        %         data = reshape(data,[frames voxels]);
        %         data = medfilt1(data,7);
        %
        %         Maxs = max(data);
        %         [~,MaxInds] = sort(-Maxs);
        %         data = data(:,MaxInds(1:min(20,size(data,2))));
        %         figure(77);
        %         subplot(311)
        %         plot(data);
        %
        %         for ind=1:size(data,2)
        %             inds = find(data(:,ind)>max(data(:,ind))/2);
        %             Stds(ind) = max(inds)-min(inds);
        %         end
        %         [~,StdInds] = sort(Stds);
        %         data = data(:,StdInds(1:min(5,size(data,2))));
        %         figure(77);
        %         subplot(312)
        %         plot(data);
        %
        %         Tails = mean(data(round(frames/5):end,:)) ./ max(data);
        %         [~,TaiInds] = sort(Tails);
        %         data = data(:,TaiInds(1:min(5,size(data,2))));
        %         figure(77);
        %         subplot(413)
        %         plot(data,'LineWidth',2);
        %         legend('1','2','3','4','5')
        
%% Zobrazení dat

% imagesc(output(:,:,1,1)-data(:,:,5,end),[0,300])
% imagesc((output(:,:,1,1)-data(:,:,5,end))./data(:,:,5,end),[0,300])
% imagesc(output(:,:,1,1)-data(:,:,5,end),[0,300])
% output111 = imrotate(output(:,:,5,1),90))
% data111 = imrotate(data(:,:,5,end),90))
% plot(squeeze(data(50,104,5,:)))
% imagesc(output(:,:,1,1)-data(:,:,5,end),[0,1])
% plot(squeeze(data(66,9,5,:)))
% plot(squeeze(data(66,97,5,:)))

%% Příprava workspace
clear all;
clc;
close all;

%% Načtení datasetů
% load('/home/karela/data/StrokeRatBrain/batches/20190610_092401_190610_STROKE_RAT_BRAIN_AMCAO9_1_1__01/experiment1/conv/data.mat');
% data_reg = load('/home/karela/data/StrokeRatBrain/batches/20190610_092401_190610_STROKE_RAT_BRAIN_AMCAO9_1_1__01/experiment1/reg/data.mat');

%load('/home/karela/data/StrokeRatBrain/batches/20190611_093556_190611_STROKE_RAT_BRAIN_AMCAO10_1_1_01/experiment1/conv/data.mat');
%data_reg = load('/home/karela/data/StrokeRatBrain/batches/20190611_093556_190611_STROKE_RAT_BRAIN_AMCAO10_1_1_01/experiment1/reg/data.mat')

%load('/home/karela/data/StrokeRatBrain/batches/20190614_094953_20190614_Stroke_RAT_BRAIN_AMCAO12_1_1_01/experiment1/conv/data.mat');
%data_reg = load('/home/karela/data/StrokeRatBrain/batches/20190614_094953_20190614_Stroke_RAT_BRAIN_AMCAO12_1_1_01/experiment1/reg/data.mat');

%load('/home/karela/data/StrokeRatBrain/batches/20190821_092549_190821_STROKE_RAT_BRAIN_AMCAO14_1_1_01/experiment1/conv/data.mat')
%data_reg = load('/home/karela/data/StrokeRatBrain/batches/20190821_092549_190821_STROKE_RAT_BRAIN_AMCAO14_1_1_01/experiment1/reg/data.mat')

% load('/home/karela/data/StrokeRatBrain/batches/20200610_101545_AMCAO_AC05_1_1_20200610_101545_AMCAO_AC05_1_1_01/experiment1/conv/data.mat');
% data_reg = load('/home/karela/data/StrokeRatBrain/batches/20200610_101545_AMCAO_AC05_1_1_20200610_101545_AMCAO_AC05_1_1_01/experiment1/reg/data.mat');

%load('/home/karela/data/StrokeRatBrain/batches/20200611_114511_AMCAO_AC08_1_1_20200611_114511_AMCAO_AC08_1_1_01/experiment1/conv/data.mat');
%data_reg = load('/home/karela/data/StrokeRatBrain/batches/20200611_114511_AMCAO_AC08_1_1_20200611_114511_AMCAO_AC08_1_1_01/experiment1/reg/data.mat');

%load('/home/karela/data/StrokeRatBrain/batches/20200713_105447_AMCAO_AC10_1_1_20200713_105447_AMCAO_AC10_1_1_01/experiment1/conv/data.mat');
%data_reg = load('/home/karela/data/StrokeRatBrain/batches/20200713_105447_AMCAO_AC10_1_1_20200713_105447_AMCAO_AC10_1_1_01/experiment1/reg/data.mat');

%% Simulovaná data
load('SyntheticEchoes_cartesian_reco');
load('SyntheticEchoes_cartesian_ReferenceMaps');

%% Angiografický dataset
% load('/nfs/bubla/data/dce-mri/upt/StrokeRatBrain/ResliceTOFonCS/20200610_101545_AMCAO_AC05_1_1/TOF_resliced/TOF_resliced.mat');
%load('/nfs/bubla/data/dce-mri/upt/StrokeRatBrain/ResliceTOFonCS/20190610_092401_190610_STROKE_RAT_BRAIN_AMCAO9_1_1/TOF_resliced/TOF_resliced.mat');
%load('/nfs/bubla/data/dce-mri/upt/StrokeRatBrain/ResliceTOFonCS/20190610_092401_190610_STROKE_RAT_BRAIN_AMCAO9_1_1/TOF_resliced/TOF_resliced.mat');

% data_reg = data_reg.data2;
% nSlices=size(data_reg,1);
% nSamples=size(data_reg,2);
% [nRows,nCols]=size(data_reg{1,1});
% data4D_reg=zeros(nRows,nCols,nSlices,nSamples);
% for sl=1:nSlices
%     for t=1:nSamples
%         data4D_reg(:,:,sl,t)=imrotate(data_reg{sl,t},90);
%     end
% end

%% Convert cell data to 4D matrix (x,y,z,time)
nSlices=size(data2,1);
nSamples=size(data2,2);
[nRows,nCols]=size(data2{1,1});
data4D=zeros(nRows,nCols,nSlices,nSamples);
for sl=1:nSlices
    for t=1:nSamples
        data4D(:,:,sl,t)=imrotate(data2{sl,t},90);
    end
end

outputFolder=[];
nSamples=size(data4D,4);
timeAxis=(0:1:nSamples-1)*info.acq.Ts/60;
Ts = info.acq.Ts;

% data4D=permute(data4D,[2 1 3 4]);

% data_reg = data4D_reg;
data = data4D;

%% Simul
data2D=DataCells2Mat(data2);
% 
mean_baseline = mean(data2D(:,:,:,1:15),4);
%% Převod pro simul
for i = 1:size(data2D,4)
% data(:,:,:,i) = (data2D(:,:,:,i)-data2D(:,:,:,1))./data2D(:,:,:,1);
data(:,:,:,i) = (data2D(:,:,:,i)-mean_baseline)./mean_baseline;
end

data4D = data;

mean_baseline = mean(SIs(:,:,1:2000),3);
%% Převod pro simul
for i = 1:size(SIs,3)
% data(:,:,:,i) = (data2D(:,:,:,i)-data2D(:,:,:,1))./data2D(:,:,:,1);
SIs_1(:,:,i) = (SIs(:,:,i)-mean_baseline)./mean_baseline;
end

%% Normalizace dat
%% Maximem
% max_data = max(data,[],'all');
% data = data./max_data;

% %AUC
% data = data./(sum(abs(data),4)/size(data,4));

% Medián
AUC_data = sum(abs(data),4)/size(data,4);
AUC_data = AUC_data(:);
AUC_data(AUC_data==0)=[];
med_data = median(AUC_data,'all');
data = data./med_data;

% data(isinf(data))=0;
% data(isnan(data))=0;

% Simul
% data = data./10;

% data = normalize(data);
% data = data./(sum(abs(data),4)/size(data,4));
% data(isinf(data))=0;
% data(isnan(data))=0;

% load('aif_test')
% [output] = fit_biexp(aif)
% figure;plot(squeeze(aif));
% hold on
% plot(output(1,1,1,8):933,bi_exp(squeeze(output),output(1,1,1,8):933))
% xlabel('Časová osa [vzorky]')
% ylabel('AIF [a.u.]')

%% Automatic detection

[output_common] = automatic_detection(data,Ts);

% rgbImage2 = cat(3, data(:,:,7,30), output_common(:,:,7), TOF_resliced(:,:,7));
% imshow(rgbImage2);

% % Visual vs angio
% 
% I1 = data_reg(:,:,8,61);
% I1 = mat2gray(I1);
% I1 = I1.*255;
% I1 = uint8(I1);
% 
% a = output_common(:,:,8);
% a = uint8(~a);
% 
% I2 = I1.*a;
% I3 = I1;
% I3(a==0)=255;
% 
% 
% I1(:,:,1) = I3;
% I1(:,:,2) = I2;
% I1(:,:,3) = I2;
% 
% figure;
% subplot 121
% imshow(I1);
% subplot 122
% TOF_resliced = imrotate(TOF_resliced,180);
% TOF_resliced = fliplr(TOF_resliced);
% imshow(TOF_resliced(:,:,8),[])

%% Simul visual

I1 = data2D(:,:,:,61);
I1 = mat2gray(I1);
I1 = I1.*255;
I1 = uint8(I1);

a = output_common(:,:,:);
im_ref = output_common(:,:,:);
a = uint8(~a);

I2 = I1.*a;
I3 = I1;
I3(a==0)=255;


I1(:,:,1) = I3;
I1(:,:,2) = I2;
I1(:,:,3) = I2;

figure;
subplot 121
imshow(I1);
subplot 122
imshow(RoisMap,[])

%% Kvant.
% Shoda
ref_map = RoisMap;
ref_map = ref_map > 37 & ref_map < 54;
ref_map = imresize(ref_map, [128 128]);

Shoda = ref_map(ref_map == im_ref);
disp(Shoda);

figure;
subplot 121
imshow(im_ref);
subplot 122
imshow(ref_map);


% Dice
dice_coeff = dice(im_ref,ref_map);
disp(dice_coeff);


% figure;
% imshowpair(output_common(:,:,6),TOF_resliced(:,:,6));

% median
% output_med5 = output_common(:,:,5);
% output_med5 = medfilt2(output_med5,[2 2]);

figure
imshow(output_common(:,:,8))

aifs_mean = zeros(size(data,3),size(data,4));

for j = 8
[rows, columns]=find(output_common(:,:,j)==1);

%% Get AIF
aif = zeros(size(data,4)*size(data,3),length(rows));
aif = [];
aifs = [];
aif_s = [];

for i = 1:length(rows)
        %aifs = squeeze(data(rows(i,1),columns(i,1),:,:));  
        aifs= squeeze(data(rows(i,1),columns(i,1),j,:));
        aif_s(i,:) = aifs; 
        %aif_s(1,1:size(data,4))=aifs(1,:);
        %for k = 1:10
            %aif_s(1,k*size(data,4)+1:(k+1)*size(data,4))=aifs(k+1,:);
        %end
        %aif(:,i) = aif_s;
end

figure;
plot(aif_s');
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')

aif_mean = mean(aif_s,1);
figure;
plot(aif_mean);
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')

% hold on
% plot(aif_ref); 
% xlabel('Časová osa [vzorky]')
% ylabel('Koncentrace k.l. [a.u.]')
aif_e = isempty(aif_mean);

if aif_e==1;
aifs_mean(j,:) = zeros(1,size(data,4)); 
else
aifs_mean(j,:) = aif_mean; 
end

% area = trapz(aif);
% aif1 = aif/area;

% for z = 1:size(aif,2)
% aif1(:,z) = aif(:,z)/max(aif(:,z));
% end
% 
% aif1(isnan(aif1))=0;

% aif2 = zeros()
% for z = 1:size(aif,2)
% aif2(:,z) = aif1(~isnan(aif1(:,z)));
% end

% figure;
% %plot(aif1*Fs);
% plot(aifs(j,:));
% 
% 
% aif_mean = mean(aifs(j,:),1);
% figure;
% %plot(aif_mean*Fs);
% plot(aif_mean);

end

aifs_mean = mean(aifs_mean,1);
figure;
plot(timeAxis,aifs_mean);
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')

%aif_09 = aifs_mean;
%aif_10_1 = aifs_mean;
%aif_12 = aifs_mean;
%aif_14 = aifs_mean;
%aif_05 = aifs_mean;
%aif_08 = aifs_mean;
% aif_10 = aifs_mean;

%% Pseudoref AIF
% % AMCAO9
% figure;
% plot(timeAxis,squeeze(data4D(34,44,8,:)));
% hold on
% plot(timeAxis,aifs_mean);
% xlabel('Časová osa [s]')
% ylabel('AIF [a.u.]')

%output_over(:,:,:) = output1(output1(:,:,5)==output2(:,:,5) & output1(:,:,5)==output3(:,:,5))

% imshowpair(rgbImage1,data(:,:,5,end));

% full_image = rgbImage1 + data(:,:,5,end);
% full_image = imrotate(full_image,90);
% imshow(full_image);

%aif1 = squeeze(data(66,97,5,:));

% for i = 1:size(output,1)
%     for j = 1:size(output,2)
%         for k = 1:size(output,3)
%         if output(i,j,k,3)>1.5
%            output1(i,j,k) = output(i,j,k,3);
%         else
%            output1(i,j,k) = 0;
%         end
%         end
%     end
% end

% output1 = output(:,:,:,3);
% output2 = zeros(size(output1,1),size(output1,2));
% 
% for i = 1:size(output1,1)
%     for j = 1:size(output1,2)
%         if output1(i,j)< 1.5
%            output2(i,j) = 0;
%         else
%            output2(i,j) = output1(i,j);
%         end
%     end
% end
% for i=1:11
%         figure;
%         imagesc(output(:,:,i,3)-data(:,:,i,end),[0,300]);
% end

%% Radovan             
        % ind = input('Which AIF number to select? ');
        % aif = data(:,best_width);
%         figure(77);
%         subplot(313)
%         plot(aif);
        
        
        %         tissue{1,1} = aif;
        %         tissue{1,2} = logical(mask);
        %         tissue{1,3} = [ Polygons{1,1}/svX*X, Polygons{1,2}/svY*Y ];
        aif={aif(:)'};
        
        info.aif.method = 'autoselect';
        info.aif.model = 'none';
        %info.aif.scaling = ScalingMethod;
        
        %% Plot AIF for next step:
        t = info.acq.Ts *(0:length(data2)-1) / 60;
        h=figure;
        size(t)
        size(aif{1})
        if length(aif{1}) < length(t)
            aif{1} = [aif{1} zeros(1,length(t)-length(aif{1}))];
        end
        if length(aif{1}) > length(t)
            aif{1} = aif{1}(1:length(t));
        end
        plot(t,aif{1}); xlabel('time [min]'); ylabel('AIF [a.u.]')
        set(h,'position',[100 100 900 400])
        %print(h,'-dpng',OutputPlotName);
        print(h,'-dsvg',OutputPlotName);
        
    case 'blind'
        %% Blind deconvolution
        
        %% Need for temporal upsampling?
        if info.acq.Ts > 1.5
            disp('Upsampling to Ts=1s')
            t = info.acq.Ts *(0:size(data2,2)-1) / 60; %min
            t_new = 0:1/60:max(t); %min (temp. sampling 1 sec)
            
            d2 = DataCells2Mat(data2);
            Zs = size(d2,3);
            Ys = size(d2,1);
            Xs = size(d2,2);
            d3 = zeros(Ys,Xs,Zs,length(t_new));
            for z=1:Zs
                for y=1:Ys
                    for x=1:Xs
                        d3(y,x,z,:) = interp1(squeeze(t),squeeze(d2(y,x,z,:)),t_new,'spline');
                    end
                end
            end
            data2 = Mat2DataCells(d3);
            info.acq.Ts = 1;
            movefile(InpFile, [InpFile(1:end-4) '_oldTs.mat']);
            save(InpFile,'data1','data2','info','-v7.3');
        end
        
        
        
        X = size(data2{1},2)
        Y = size(data2{1},1)
        Z = size(data2,1)
        
        %% Create masks and extract tissue curves to 'tissue' - for voxels
        tissue = {};
        frames = size(data2,2);
        
        %         %% Clustering
        %         mask = zeros(Y,X,Z);
        %         NumPoly = size(Polygons,1);
        %         % union of masks
        %         for ind=1:NumPoly
        %             mask(:,:,Polygons{ind,5}) = mask(:,:,Polygons{ind,5}) + poly2mask(Polygons{ind,1}/svX*X,Polygons{ind,2}/svY*Y,size(data2{1},1),size(data2{1},2));
        %         end
        %         mask(mask>0)=1;
        %
        %         curves = zeros(sum(sum(sum(mask))),frames);
        %         for f=1:frames
        %             cube = DataCells2Mat(data2(:,f));
        %             curves(:,f) = cube(mask>0);
        %         end
        %
        %         NumOfChannels = 10;
        %         [~, curves2]=kclusterSim(curves,NumOfChannels);
        %         NumOfChannels = size(curves2,1);
        %
        %         %% extract clusters with zero curve
        %         curves3=[];
        %         for ii=1:NumOfChannels
        %             if sum(curves2(ii,:))>0
        %                 curves3 = [curves3; curves2(ii,:)];
        %             end
        %         end
        %         NumOfChannels = size(curves2,1);
        %
        %         for curve=1:NumOfChannels
        %             tissue{curve,1} = curves3(curve,:);
        %             tissue{curve,2} = logical(mask);
        %             tissue{curve,3} = [ Polygons{1,1}/svX*X, Polygons{1,2}/svY*Y ];
        %             figure(1010); plot(tissue{curve,1}); hold on
        %         end
        %         hold off
        %         %pause
        
        if strcmp(BlindModel,'TH+TV')
            
            %% single manually entered ROIs - extract each voxel separately
            NumPoly = size(Polygons,1);
            Polygons
            for ind=1:1 %NumPoly  %%% only the first ROI !!!!!!!!!!!!!
                %mask = poly2mask(Polygons{ind,1}/zvetseni,Polygons{ind,2}/zvetseni,size(data2{1},1),size(data2{1},2));
                ind
                Polygons{ind,1}
                Polygons{ind,2}
                %% FOR ROIs
                %mask = poly2mask(Polygons{ind,1}/svX*X,Polygons{ind,2}/svY*Y, Y,X);
                mask = zeros(Y,X,Z);
                mask(:,:,Polygons{ind,5}) = poly2mask(Polygons{ind,1}/svX*X,Polygons{ind,2}/svY*Y, Y,X);
                [row,col] = find(mask(:,:,Polygons{ind,5}));
                for poin = 1:length(row)
                    maskP = zeros(Y,X);
                    maskP(row(poin),col(poin)) = 1;
                    tissue{poin,2} = logical(maskP);
                    tissue{poin,3} = [ Polygons{ind,1}/svX*X, Polygons{ind,2}/svY*Y ];
                    curve = zeros(frames,1);
                    for f=1:frames
                        curve(f) = data2{Polygons{ind,5},f}(row(poin),col(poin));
                    end
                    tissue{poin,1} = curve;
                end
            end
            
        else
            
            %% single manually entered ROIs
            NumPoly = size(Polygons,1);
            NotNanInd = 1;
            for ind=1:NumPoly
                %mask = poly2mask(Polygons{ind,1}/zvetseni,Polygons{ind,2}/zvetseni,size(data2{1},1),size(data2{1},2));
                ind
                Polygons{ind,1}
                Polygons{ind,2}
                %% FOR ROIs
                %mask = poly2mask(Polygons{ind,1}/svX*X,Polygons{ind,2}/svY*Y, Y,X);
                mask = zeros(Y,X,Z);
                mask(:,:,Polygons{ind,5}) = poly2mask(Polygons{ind,1}/svX*X,Polygons{ind,2}/svY*Y, Y,X);
                %% FOR POINTs
                %             mask = zeros(size(data2{1}));
                %             mask(round(Polygons{ind,2}(1)/svY*Y), round(Polygons{ind,1}(1)/svX*X))  = 1;
                curve = zeros(frames,1);
                for f=1:frames
                    curve(f) = sum(sum(sum(mask .* DataCells2Mat(data2(:,f)) ))) / sum(sum(sum(mask)));
                end
                if isfinite(sum(curve)) % ROI not empty:
                    tissue{NotNanInd,1} = curve;
                    tissue{NotNanInd,2} = logical(mask);
                    tissue{NotNanInd,3} = [ Polygons{ind,1}/svX*X, Polygons{ind,2}/svY*Y ];
                    NotNanInd = NotNanInd + 1;
                end
            end
            
        end
        
        
        %% parameters of blind deconvolution:
        
        info.codewords{1,size(info.codewords,2)+1}='aif';
        aif = cell(1,1);
        %selected_model =  'Gamma-variate Tri-exponencial';
        %selected_model =  'Parker';
        %selected_model =  'DiBellaLim';
        %selected_model =  'DiBellaLim2';
        selected_model =  AIFModel;
%        t = info.acq.Ts *(0:length(tissue{1,1})-1) / 60;
        
        AIFinits = [...
            0.5174    0.4493    2.1692    3.1172    0.0703    1.5942    0.0822; % Mean blindly estimated Magnevist AIF
            %1.1608    0.0564    0.1502    5.0731    0.1547    0.9921    0.4729; % reference
            %0.1230    0.1808    1.7559    0.0840    5.4305    1.8581    0.4000; % mine
            %1.4742    0.8100         0    2.1400    0.0900    1.1150    0.1000; % Ragan
            %1.2000    0.4700         0    1.2000    0.0500    0.6250    0.1000; % Chockley
            %0.9200    0.4200         0    2.3000    0.0500    1.1750    0.1000; % Heilmann
            ];
        crit = zeros(1,size(AIFinits,1));
        
        
        %% call blind deconvolution:
        
        for InitAifInd = 1:size(AIFinits,1)
            AIFinit = AIFinits(InitAifInd,:);
            %[aif{InitAifInd} info.aifpar crit(InitAifInd)]= multich_DCATH_alt_batch_moreiter_muscle(selected_model,tissue,t, 0 ,AIFinit,PBfile);
            %[aif{InitAifInd} info.aifpar crit(InitAifInd) tissue]= multich_DCATH_alt_batch2(selected_model,tissue,t, 0 ,AIFinit,PBfile);
            %[aif{InitAifInd} info.aifpar crit(InitAifInd) tissue]= multich_DCATH_alt_batch3_BBB(selected_model,tissue,t, 0 ,AIFinit,NumBBB,BlindModel,PBfile);
            %[aif{InitAifInd} info.aifpar crit(InitAifInd) tissue]= multich_DCATH_alt_batch3_BBB02(selected_model,tissue,t, 0 ,AIFinit,NumBBB,BlindModel,PBfile);  % nested, not well functionaing
            %[aif{InitAifInd} info.aifpar crit(InitAifInd) tissue]= multich_DCATH_alt_batch3_BBB03(selected_model,tissue,t, 0,AIFinit,NumBBB,BlindModel,PBfile);  % alternating, not well functionaing
            %[aif{InitAifInd} info.aifpar crit(InitAifInd) tissue]= multich_DCATH_alt_batch3_BBB04(selected_model,tissue,t, NumBBB,BlindModel,PBfile);  % alternating, using nonblind dec. of Michal
            %[aif{InitAifInd} info.aifpar crit(InitAifInd) tissue]= multich_DCATH_alt_batch3_BBB05(selected_model,tissue,t, NumBBB,BlindModel,info,PBfile);  % alternating, using nonblind dec. of Michal but constraining MTT (a new model called 'dcs' used in Michal's code)
            
            %             if strcmp(BlindModel,'TH+TV')
            %% regularized
%            [aif{InitAifInd} info.aifpar crit(InitAifInd) tissue]= multich_DCATH_alt_batch3_BBB05reg(selected_model,tissue,t, NumBBB,BlindModel, info, PBfile);  % alternating, using regularized nonblind dec. of Michal
            %             else
            %                 %% standard
            %                 [aif{InitAifInd} info.aifpar crit(InitAifInd) tissue]= multich_DCATH_alt_batch3_BBB05(selected_model,tissue,t, NumBBB,BlindModel, info, PBfile);  % alternating, using nonblind dec. of Michal but constraining MTT (a new model called 'dcs' used in Michal's code)
            %             end
        end
        crit
        [m, optind] = min(crit);
        aif{1} = aif{optind};
        aif = aif(1);
        
        info.aif.method = 'BlindDec'
        info.aif.model = selected_model;
        %info.aif.scaling = ScalingMethod;
        
        %% Plot AIF for next step:
        h=figure;
        
        subplot(1,2,1);
%        plot(t,aif{1}); xlabel('time [min]'); ylabel('AIF [a.u.]')
        
        subplot(1,2,2);
        Colors = jet(size(tissue,1));
        hold on;
        for ii=1:size(tissue,1)
            plot(t,tissue{ii,1},':.','Color',Colors(ii,:));
            plot(t,tissue{ii,end}/sum(tissue{ii,end})*sum(tissue{ii,1}),'-','Color',Colors(ii,:));
        end
        hold off
        
        set(h,'position',[100 100 900 400]);
        %print(h,'-dpng',OutputPlotName);
        print(h,'-dsvg',OutputPlotName);
        
        
        if strcmp(BlindModel,'TH+TV')
            [maps names units ranges] = GetParamMapsGeneral_TH(tissue);
            save(MapFile, 'maps', 'names', 'units', 'ranges');
            for zz = 1:Z
                for map=1:length(maps)
                    nm = names{map};
                    nm = nm(isletter(nm));
                    imwrite(imresize(  ((maps{map}(:,:,zz))-ranges{map}(1)) / (ranges{map}(2)-ranges{map}(1))  *64,  [svY svX]),...
                        jet, fullfile(ViewerParDir, [sprintf('Slice%02d_',zz) nm '.jpg']),'Quality',100);
                end
            end
        end
        
        
end

%% save
save(OutFile,'info','aif','tissue');

PB(PBfile,100);


