function T1_lsf = analyzeT1(data1,info,SeqType,rectROI,setFrame)
%% MAIN FUNCTION
% Function T1_lsf calculates the T1 map from MRI data set.
% The input parameters are:
%   o data1 - data set of MRI images
%   o info - structure that contains information about aquizition
%   o SeqType - information about the sequence type
%             - 1 for FISP FID, 2 for Intragate FLASH
%   o rectROI - binary mask of rectangular ROI
%   o setFrame - vector [q,s], where s is the chosen slice and q is the
%                chosen cardiac slice
% The output parameter
%   o T1_lsf - matrix of calculated T1 values

% preset of some variables
siz_x = size(data1);
siz_y = size(data1{1,1}{1,1});
siz_z = size(data1{1});
T1_lsf = zeros(siz_y(1),siz_y(2),siz_z(2));
M0 = zeros(siz_y(1),siz_y(2),siz_z(2));
vec = zeros(1,siz_x(2));
TI = zeros(1,siz_x(2));
FA = zeros(1,siz_x(2));
s = setFrame(1);
q = setFrame(2);
if SeqType == 1 % setting the TI or FA vector
    for l = 1:siz_x(2)
        TI(l) = info.acq.TI(l);
    end
else
    for l = 1:siz_x(2)
        FA(l) = info.acq.FA(l);
    end
end
for i = 1:siz_y(1)
    for j = 1:siz_y(2)
        if rectROI(i,j)
            if SeqType == 1 % for FISP FID
                for k = 1:siz_x(2)
                    TI(k) = info.acq.TI(k);
                    vec(1,k) = data1{1,k}{q,s}(i,j);
                end
                T10 = [1.7e6 1500]; % initial values
                % using lsfcurvefit function for fitting
                [est] = lsqcurvefit(@TI_fun,T10,TI,vec,[0 0],[1e10 2500]); 
                T1_lsf(i,j) = est(2);
                M0(i,j) = est(1);
            else % for Intragate FLASH
                for k = 1:siz_x(2)
                    FA(k) = info.acq.FA(k);
                    vec(1,k) = data1{1,k}{q,s}(i,j);
                end
                FA = FA*(pi/180); % conversion from ° to rad
                T10 = [5e6 500]; % initial values
                % using lsfcurvefit function for fitting
                est = lsqcurvefit(@FA_fun,T10,FA,vec,[0 0],[1e10 2500]); 
                T1_lsf(i,j) = est(2);
                M0(i,j) = est(1);
            end
        end
    end
end
% deleting rows and columns with all zeros
T1_lsf(:,all(T1_lsf==0,1))=[]; 
T1_lsf(all(T1_lsf==0,2),:)=[];
end

%% ADDITIONAL FUNCTIONS
function f = TI_fun(x,TI)
% Function TI_fun evaluates the FID value of IR sequence.
% The input parameters
%   o x - vector [M0 T1]
%   o TI - vector of TI values
% The output parameters
%   o f - vector of FID values
f = abs(x(1)*(1-2*exp(-TI/x(2))));
end

function f = FA_fun(x,alfa)
% Function TI_fun evaluates the FID value by Ernst equation.
% The input parameters
%   o x - vector [M0 T1]
%   o alfa - vector of FA values in rad
% The output parameters
%   o f - vector of FID values
global info
% retrieving the repetition time from info structure
TR = info.acq.TR;
M = x(1);
E1 = exp(-TR/x(2));
f = abs((M*sin(alfa))*(1-E1)./(1-E1*cos(alfa)));
end