function [AUCfiltered,ROUfiltered,TF_ROU] = AUCandROU(data4D,info,outputFolder,AUCtreshold,ROUtreshold,ROUevalFromSample,visu,saveFile)
%thresholds
% AUCtreshold=0.97;%10 percent nejvetsich ponechat dle clanku https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.20759
% ROUtreshold=0.75; %procento nejmene roztrepanych krivek dle clanku https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.20759
% ROUevalFromSample=30; %calculate roughness from 30th sample until the end of CTC sequence

%% create and display AVG matrix
data4DAVG=mean(data4D,4);
minAVG=min(data4DAVG(:));
maxAVG=max(data4DAVG(:));

%% calculate AUC
AUCfull=sum(data4D,4);
minAUC=min(AUCfull(:));
maxAUC=max(AUCfull(:));

nVzorku=numel(AUCfull);
AUC1dsorted=sort(AUCfull(:),'ascend');
treshIdx=ceil(nVzorku*AUCtreshold);
threshValue=AUC1dsorted(treshIdx);
TF_AUC=false(size(AUCfull));
TF_AUC(AUCfull>threshValue)=true;
AUCfiltered=AUCfull.*TF_AUC;
% disp(['propusteno procent voxelu ' num2str((sum(TF_AUC(:))/numel(TF_AUC))*100) ]);
% num2str((sum(TF_AUC(:))/numel(TF_AUC))*100
% disp(['Reduction factor based on AUC ' num2str(treshIdx/nVzorku*100) char(37)]);

% figure('name','AUC calculated');
% montage(AUCfull,'Size',[4 4],'DisplayRange',[0 maxAUC*0.5]);
figure('name','AUC mask');
montage(TF_AUC,'Size',[4 4],'DisplayRange',[0 1]);


%% Calculate roughness proper
% close all force
% normalize data to CTC have unity AUC
data4DNorm=data4D./repmat(sum(data4D,4),[1 1 1  size(data4D,4)]);
dif=2; %2nd difference
procDim=4; %across time domain
data4D2nddiff = diff(data4DNorm,dif,procDim); %druha derivace pres casovou 4. dimenzi
ROUfull   = sum(data4D2nddiff(:,:,:,ROUevalFromSample:end).^2,procDim);%normalizovana 3D matice ROUGHNESS
ROUfull(isnan(ROUfull))=inf; %remove NaNs

nVzorkuPoRedukciAUC=sum(TF_AUC(:));
treshIdx=ceil(nVzorkuPoRedukciAUC*ROUtreshold); %10 percent


% disp(['Reduction factor based on ROU ' num2str(treshIdx/nVzorkuPoRedukciAUC*100) char(37)]);
% disp(['Complete reduction factor based on AUC and ROU ' num2str(treshIdx/numel(AUCfull)*100) char(37)]);

ROUreduced=ROUfull(TF_AUC);
ROU1drecucedSorted=sort(ROUreduced(:),'ascend');
threshValue=ROU1drecucedSorted(treshIdx);

TF_ROU=false(size(ROUfull));
TF_ROU(ROUfull<threshValue)=true;

TF_ROU=TF_ROU & TF_AUC;

ROUfiltered=ROUfull.*TF_ROU;

if visu
    figure('name','AUC filtered');
    montage(AUCfiltered,'Size',[4 4],'DisplayRange',[0 max(AUCfiltered(:))]);
    
    % figure('name','Roughness calculated');
    % montage(ROUfull,'Size', [4 4],'DisplayRange', [0 1]);
    
    figure('name','Roughness mask');
    montage(TF_ROU,'Size', [4 4],'DisplayRange', [0 1]);
    figure('name','Roughness filtered');
    montage(ROUfiltered,'Size', [4 4],'DisplayRange', [0 max(ROUfiltered(:))]);
    
end
if saveFile
    disp('Saving roughness matrix');
    save(fullfile(outputFolder,'ROUmatrix.mat'),...
        'AUCfiltered','ROUfiltered','TF_ROU','data4D','info','AUCtreshold','ROUtreshold','ROUevalFromSample');
end
disp('AUC a ROU filtrace dokončena');
end

