function [QALAS, profileVec] = fitROI2step(QALAS,fitMode,doProfile,saveFolder)
% fitROI is a function for estimating R1 R2 and Alpha for a Vector of values
% for estimaing a enterir volume use fitQALAS
% 
%  [QALAS, R1roi, R2roi, c, Alpharoi ] = fitROI(QALAS,doProfile)
%
% Writen by Anders Tisell, 2024-12-11

testR1_1 = single([0.01 0.05 0.1 0.25 0.5:0.5:2 2.25 2.5 3 4 6]);
testR2_1 = testR1_1*10;
testAlpha_1 = single(1);
testBeta_1 = single(1);%single(0.85);%

lengthR1_1 = length(testR1_1);
lengthR2_1 = length(testR2_1);
lengthAlpha_1 = length(testAlpha_1);
lengthBeta_1 = length(testBeta_1);

R1vec_1 = repmat(testR1_1',lengthR2_1*lengthAlpha_1*lengthBeta_1,1);
R2vec_1 = repmat(sort(repmat(testR2_1',lengthR1_1,1)),lengthAlpha_1*lengthBeta_1,1);
Alphavec_1 = repmat(sort(repmat(testAlpha_1',lengthR1_1*lengthR2_1,1)),lengthBeta_1,1);
Betavec_1 = sort(repmat(testBeta_1',lengthR2_1*lengthR1_1*lengthAlpha_1,1));

testR1_add = single(0.3*randn(1,60));
testR2_add = single(3*randn(1,60));
testAlpha_2 = single(1);%single(0.7:0.01:1.3); %ssingle(0.7:0.01:1.3);%
testBeta_2 = single(1);%single(0.85); %single(0.7:0.01:1.3);%single(0.85);%single(0.7:0.01:1.3);

lengthR1_add = length(testR1_add);
lengthR2_add = length(testR2_add);
lengthAlpha_2 = length(testAlpha_2);
lengthBeta_2 = length(testBeta_2);

R1vec_2 = repmat(testR1_add',lengthR2_add*lengthAlpha_2*lengthBeta_2,1);
R2vec_2 = repmat(sort(repmat(testR2_add',lengthR1_add,1)),lengthAlpha_2*lengthBeta_2,1);
Alphavec_2 = repmat(sort(repmat(testAlpha_2',lengthR1_add*lengthR2_add,1)),lengthBeta_2,1);
Betavec_2 = sort(repmat(testBeta_2',lengthR2_add*lengthR1_add*lengthAlpha_2,1));

for ROIidx = 1:length(QALAS.ROIvalue)
     

    disp(['Rough Fit ' QALAS.ROIname{ROIidx}])

    disp('Init fit data')
    ROI_N = length(QALAS.ROIvalue{ROIidx});
    %   if isempty(doProfile) || doProfile == 0
    R1roi_1 = shiftdim(repmat(R1vec_1,1,ROI_N),1);
    R2roi_1 = shiftdim(repmat(R2vec_1,1,ROI_N),1);
    Alpharoi_1 = shiftdim(repmat(Alphavec_1,1,ROI_N ),1);
    Betaroi_1 = shiftdim(repmat(Betavec_1,1,ROI_N),1);

    tic
    disp('Rough fit QALAS')

    %% Fit with Beta factor

    estQALAS = squeeze(modelQALAS(R1roi_1,R2roi_1,QALAS.info,Alpharoi_1,Betaroi_1));
    QALASroi = single(repmat(QALAS.ROIvalue{ROIidx},1,1,size(R1roi_1,2)));
    QALASroi = permute(QALASroi,[1 3 2]);
    QALASroiMean = repmat(mean(QALASroi,1),ROI_N,1,1);
    c_1 = mean(QALASroiMean./estQALAS,3);

    disp('Calculate min')

    minC = 50;
    maxC = 2000;
    c_1(c_1<50) = minC;
    c_1(c_1>2000) = maxC;
   % testM0_1 = minC:maxC;
    %

    meanROIresidual = sum(sqrt((c_1.*estQALAS - QALASroiMean).^2),3);
    [~, minIDX] =  min(meanROIresidual,[],2,'linear');

    residual = sum(sqrt((c_1.*estQALAS - QALASroi).^2),3);
    
    estM0_1 = c_1(minIDX);
    estR1_1 = R1roi_1(minIDX);
    estR2_1 = R2roi_1(minIDX);
    estAlpha_1 = Alpharoi_1(minIDX);   
    estBeta_1 = Betaroi_1(minIDX);

%% Fin tune the fitting 

    R1roi_2 = shiftdim(repmat(R1vec_2,1,ROI_N),1) + shiftdim(repmat(estR1_1,1,length(R1vec_2),1));
    R2roi_2 = shiftdim(repmat(R2vec_2,1,ROI_N),1) + shiftdim(repmat(estR2_1,1,length(R2vec_2),1));
    Alpharoi_2 = shiftdim(repmat(Alphavec_2,1,ROI_N ),1);% + shiftdim(repmat(estAlpha_1,1,length(Alphavec_2),1));
    Betaroi_2 = shiftdim(repmat(Betavec_2,1,ROI_N),1);% + shiftdim(repmat(estBeta_1,1,length(Betavec_2),1));
    
    R1roi_2(R1roi_2 < 0) = 0;
    R2roi_2(R2roi_2<0) = 0;

    if fitMode == 1

        estQALAS = squeeze(modelQALAS(R1roi_2,R2roi_2,QALAS.info,Alpharoi_2,Betaroi_2));

  % Get all data
        QALASroi = single(repmat(QALAS.ROIvalue{ROIidx},1,1,size(R1roi_2,2)));
        QALASroi = permute(QALASroi,[1 3 2]);

        c = mean(QALASroi./estQALAS,3);
        
        disp('Calculate min')

        minC = 50;
        maxC = 2000;
        c(c<50) = minC;
        c(c>2000) = maxC;
        testM0_2 = minC:maxC;

        residual = sum(sqrt((c.*estQALAS - QALASroi).^2),3);
        [~, minIDX] = min(residual,[],2,'linear');
        meanROIresidual = [];
    elseif fitMode == 2

        estQALAS = squeeze(modelQALAS(R1roi_2,R2roi_2,QALAS.info,Alpharoi_2,Betaroi_2));

        QALASroi = single(repmat(QALAS.ROIvalue{ROIidx},1,1,size(R1roi_2,2)));
        QALASroi = permute(QALASroi,[1 3 2]);
        c = mean(QALASroi./estQALAS,3);
        %QALASroiMean = repmat(mean(QALASroi,1),ROI_N,1,1);
        %c = mean(QALASroiMean./estQALAS,3);

        disp('Calculate min assuming one R1, R2 and B1 for the ROI but variable c')

        minC = 50;
        maxC = 2000;
        c(c<50) = minC;
        c(c>2000) = maxC;
        testM0_2 = minC:maxC;
        %
        %meanROIresidual = sum(sqrt((c.*estQALAS - QALASroiMean).^2),3);
        %[~, minIDX] =  min(meanROIresidual,[],2,'linear');

        residual = sum(sqrt((c.*estQALAS - QALASroi).^2),3);
        [ResidualRows, residualCols] = size(residual);
        FitError = repmat(sum(residual,1),[ResidualRows, 1]);


        meanROIresidual = mean(residual,1);
        [~, minIDX] = min(FitError,[],2,'linear');
    elseif fitMode == 3

        estQALAS = squeeze(modelQALAS(R1roi_2,R2roi_2,QALAS.info,Alpharoi_2,Betaroi_2));

        QALASroi = single(repmat(QALAS.ROIvalue{ROIidx},1,1,size(R1roi_2,2)));
        QALASroi = permute(QALASroi,[1 3 2]);
        
        QALASroiMean = repmat(mean(QALASroi,1),ROI_N,1,1);
        c = mean(QALASroiMean./estQALAS,3);

        disp('Calculate min avrege ROI signal')

        minC = 25;
        maxC = 2000;
        c(c<50) = minC;
        c(c>2000) = maxC;
        testM0_2 = minC:maxC;
% 
        meanROIresidual = sum(sqrt((c.*estQALAS - QALASroiMean).^2),3);
        [~, minIDX] =  min(meanROIresidual,[],2,'linear');
        residual = sum(sqrt((c.*estQALAS - QALASroi).^2),3);

    end

    toc
    tic
    disp('Write QALAS struct' )
    estR1 = R1roi_2(minIDX);
    estR2 = R2roi_2(minIDX);
    estAlpha = Alpharoi_2(minIDX);
    estM0 = c(minIDX);
    estBeta = Betaroi_2(minIDX);

    estError = residual(minIDX);

    QALAS.estROI{ROIidx}.R1 = estR1;
    QALAS.estROI{ROIidx}.R2 = estR2;
    QALAS.estROI{ROIidx}.Alpha = estAlpha;
    QALAS.estROI{ROIidx}.M0 = estM0;
    QALAS.estROI{ROIidx}.Beta = estBeta;
    QALAS.estROI{ROIidx}.Error = estError;
   

    BestFitQALAS = squeeze(modelQALAS(estR1,estR2,QALAS.info,estAlpha,estBeta));
    QALASsignalInROI = BestFitQALAS.*repmat(estM0,1,5);

    QALAS.estROIsignal{ROIidx} = QALASsignalInROI;
    QALAS.estROIresidual{ROIidx} = QALAS.ROIvalue{ROIidx} - QALASsignalInROI;

    toc
    %% Plot data
    tic
    disp('Plot data')
    if fitMode == 1
        testR1 = linspace(mean(estR1)-std(estR1)*4,mean(estR1)+std(estR1)*4,lengthR1_add);
        testR2 = linspace(mean(estR2)-std(estR2)*4,mean(estR2)+std(estR2)*4,lengthR2_add);
    elseif fitMode == 3 || fitMode == 4

        testR1 = mean(estR1)-0.2:0.01:mean(estR1)+0.2;
        testR2 = mean(estR2)-2:0.1:mean(estR2)+2;
    end

    plotQALASfit(QALAS,ROIidx,testR1,testR2,testAlpha_2,testM0_2,testBeta_2,fitMode,saveFolder)

    toc

    %% Profiles
    if doProfile
        disp('Generate Profiles')
        tic
        %profileVec = profilelikelihoodQALAS(QALAS,ROIidx,ROI_N,residual,R1roi,R2roi,Alpharoi,Betaroi,c,testR1,testR2,testAlpha,testM0,testBeta,fitMode,saveFolder);
        if isempty(meanROIresidual) 
            meanROIresidual = zeros(1,size(residual,2));
            profileVec = profilelikelihoodQALAS(QALAS,ROIidx,ROI_N,residual,meanROIresidual,R1roi_2,R2roi_2,Alpharoi_2,Betaroi_2,c,testR1,testR2,testAlpha_2,testM0_2,testBeta_2,fitMode,saveFolder);
         %   correlationsQALAS(QALAS,ROIidx,ROI_N,residual,meanROIresidual,R1roi_2,R2roi_2,Alpharoi_2,Betaroi_2,c,fitMode,saveFolder)
        else
            profileVec = profilelikelihoodQALAS(QALAS,ROIidx,ROI_N,residual,meanROIresidual(1,:),R1roi_2,R2roi_2,Alpharoi_2,Betaroi_2,c,testR1,testR2,testAlpha_2,testM0_2,testBeta_2,fitMode,saveFolder);
        end
        toc
    end
end
if isfolder(saveFolder)
    %save([saveFolder 'QALAS.mat'],'QALAS')
    save([saveFolder 'QALAS_' QALAS.patientCode '_' QALAS.examName '_' QALAS.GD '_' QALAS.protocol '.mat'],'QALAS')

end

end