function [data2, info] = SI2con_SRFLASH_whole(data1, data2, info, maska, skip, PhaseEncSteps, TR_real, precontrast, doT1map, doDynSeq)
% test conversion to a value linear w.r.t. concentration on a ROI
% input args:
%   data1, data2, info - input data
%   bw - mask of the testing roi
%   skip - how many frames to skip
%   PhaseEncSteps - number of phase encoding steps before the acquisition of the middle
%     k-space line
%   TR_real - TR of the acquisition sequence in [ms]
%   precontrast - number of precontrast frames to subtract baseline (includes the frames to be skipped)




%% precontrast sequences - extract imagesc for estimation of M0 and E10

seqs = info.acq.kM0Range(2)-info.acq.kM0Range(1)+1;   % number of precontrast seqs.
slices = info.sections;

tis = 1:3000; % TI axis

%% get image size, and number of frames
frames = 0;
for s=1:seqs
    f = length(data1{info.acq.kM0Range(1)+s-1});
    if f>frames
        frames=f;
    end
end
X = size(data1{1}{1},2);
Y = size(data1{1}{1},1);


%%
numframes = 0; % number of all precontrast images
xdata = zeros(0);
for seq = info.acq.kM0Range(1):info.acq.kM0Range(2) %1:seqs
    frn = length(data1{seq}) - skip; % number of frames of the given sequence - skipframes
    xdata = [xdata double(info.acq.TI(seq)*ones(1,frn)) ];
    numframes = numframes + frn;
end

img = zeros(Y,X,numframes);   % sequence , frame
% SI = zeros(numframes);
ind = 1;
FA = info.acq.FA(info.acq.kM0Range(1));
for seq = info.acq.kM0Range(1):info.acq.kM0Range(2) %1:seqs
    frn = length(data1{seq}) - skip; % number of frames of the given sequence - skipframes
    for frame=1:frn
        img(:,:,ind) = double( data1{seq}{frame+skip} );
        ind = ind + 1;
    end
    if info.acq.FA(seq) ~= FA
        disp('Error: different FA values in precontrast acquisition! ')
        return
    end
end


%% curve fitting to estimate E10 and M0 for complete image:

if doT1map
    % E10m = zeros(Y,X);
    M0m = zeros(Y,X);
    T10m = zeros(Y,X);
    nm1 = PhaseEncSteps;    % n-1
    warning off
    for y=1:Y %40:80 %1:Y %50:150
        for x=1:X %40:80 %1:X %50:150
            if maska(y,x)
                ydata_c = double(img(y,x,:));
                ydata = shiftdim(ydata_c,2)';
                
                M0_odhad = max(ydata);
                
                %[T10 M0] = SI2con_fit_SRFLASH(xdata, ydata, TR_real, FA, nm1, M0_odhad);
                [T10 M0] = SI2con_fit_SRFLASH(xdata, ydata, TR_real, FA, nm1);
                
                T10m(y,x)=T10;
                M0m(y,x)=M0;
                
                E10 = exp(-tis/T10);
                TR = TR_real;
                E = exp(-TR/T10);
                a = cos(FA/180*pi) .* E;
                
                yy = M0 * ( (1-E10).*(a.^nm1) + (1-E).*(1-a.^nm1)./(1-a));
                
                %figure(1);plot(xdata, ydata, 'o', tis,y,'-');
                figure(100); hold off
                plot(xdata, ydata, 'o', tis,yy);
                title(sprintf('M0:%f,  T10:%f, x:%d, y:%d',M0, T10, x, y))
            else
                T10m(y,x)=0;
                M0m(y,x)=0;
            end
        end
    end
    warning on
    
    % save(sprintf(T10mapName, input), 'T10m');
    % %save(sprintf(E10mapName, slice), 'E10m');
    % save(sprintf(M0mapName, input), 'M0m');
    
    info.maps.T10 = T10m;
    info.maps.M0  = M0m;
    info.maps.T10roi = maska;
end

figure(100); imagesc(info.maps.T10); colorbar; title('T10')
figure(101); imagesc(info.maps.M0); colorbar; title('M0')




%% Signal conversion of the dynamic sequence:

if doDynSeq
    
    %parameters
    
    TR = TR_real;
    FA = info.acq.FA(end);
    TI = info.acq.TI(end);
    nm1 = 0;    % n-1
    dyn_frames = length(data2);
    maska = logical(maska);
    
    if ~nm1
        try matlabpool open; end
    end

    parts = 1;
    part_frames = ceil(dyn_frames/parts);
    warning off
    for part = 1:parts
        parfor frame = (part_frames*(part-1)+1) : min((part_frames*part),dyn_frames)
        %for frame = 50:51
        %parfor frame = 1:dyn_frames
            
            im = double(data2{frame});
            figure(6);imagesc(im);colormap(gray);

            S = im(maska);
            M0 = info.maps.M0(maska);

            if ~nm1
                R1 = -1/TI * log(1-S./M0);
                im(maska) = R1;
            else
                T1_weight = 1000000; %100000000;
                par0 = T1_weight*exp(-1./ (ones(size(S))*300) );
                
                options=optimset(...
                    'Display','iter',...
                    'UseParallel','never',...
                    'DerivativeCheck','off','GradObj','on');
                
                myfun = @(par)final_fit(par,S, TR, TI, FA, nm1, M0, T1_weight);
                %[par,fval,exitflag,output,lambda,grad] = fmincon(myfun,par0,[],[],[],[],[zeros(size(S))],[T1_weight*ones(size(S))],[],options);
                %par = fmincon(myfun,par0,[],[],[],[],[zeros(size(S))],[T1_weight*ones(size(S))],[],options);
                par = fminimax(myfun,par0,[],[],[],[],[zeros(size(S))],[T1_weight*ones(size(S))],[],options);
                
                T1 = -1./log(par./T1_weight);
                im(maska) = 1./T1;
            end
            
            
            im=real(im);
            im(im<0)=0;
            im(~isfinite(im))=0;
            im(~maska) = 0
            data2{frame} = im;
            
            disp(sprintf('frame: %d / %d',frame,dyn_frames));
            figure(7);imagesc(data2{frame},[0.0005 0.01]);colormap(gray);title(sprintf('R1, frame: %d / %d',frame,dyn_frames))
            %figure(8);imagesc(1./data2{frame},[100 2000]);colormap(gray);title('T1')

        end
        
    end
    warning on
            
    if ~nm1
        try matlabpool close; end
    end
     
    %% subtract base line:
    baseline = zeros(size(data2{1}));
    for frame = skip+1 : precontrast
        baseline = baseline + double(data2{frame});
    end
    baseline = baseline / (precontrast-skip);
    
    for frame = 1 : length(data2)
        data2{frame} = single( data2{frame} - baseline );
    end
    
    
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ssd, g] = final_fit(par, S, TR, TI, FA, nm1, M0, T1_weight)
pom = 1;

EE = par./T1_weight;  % exp(-1/T10)

E = EE.^(TR/pom); %exp(-TR/T10);
a = cos(FA/180*pi) .* E;
E10 = EE.^(TI/pom); %exp(-TI/T10);

term1 = double(  (1-E10).*(a.^nm1) + (1-E).*(1-a.^nm1)./(1-a)  );
y = M0 .* term1;

ssd = (y-S).^2;
% ssd = mean2((y-S).^2);
% if ~isfinite(ssd)
%     ssd = 1e6;
% end

%% Derivatives:
p = par;
if nargout > 1
    % dy/dp:
    der = -M0.*((TI.*(p./T1_weight).^(TI - 1).*(cos((pi.*FA)./180).*(p./T1_weight).^TR).^nm1)./T1_weight + (TR.*((cos((pi.*FA)./180).*(p./T1_weight).^TR).^nm1 - 1).*(p./T1_weight).^(TR - 1))./(T1_weight.*(cos((pi.*FA)./180).*(p./T1_weight).^TR - 1)) + (TR.*nm1.*cos((pi.*FA)./180).*(p./T1_weight).^(TR - 1).*((p./T1_weight).^TI - 1).*(cos((pi.*FA)./180).*(p./T1_weight).^TR).^(nm1 - 1))./T1_weight - (TR.*cos((pi.*FA)./180).*((cos((pi.*FA)./180).*(p./T1_weight).^TR).^nm1 - 1).*(p./T1_weight).^(TR - 1).*((p./T1_weight).^TR - 1))./(T1_weight.*(cos((pi.*FA)./180).*(p./T1_weight).^TR - 1).^2) + (TR.*nm1.*cos((pi.*FA)./180).*(p./T1_weight).^(TR - 1).*((p./T1_weight).^TR - 1).*(cos((pi.*FA)./180).*(p./T1_weight).^TR).^(nm1 - 1))./(T1_weight.*(cos((pi.*FA)./180).*(p./T1_weight).^TR - 1)));
    g = diag(2*(y-S).* der);
    % %    if isfinite(der)
    %        g = 2*(y-S).* der;
    % %    else
    % %        g = 2*(y-S).* 1e6;
    % %    end
    %    g = 2*(y-S).* der;
    %
    % %    if ~isfinite(sum(sum(g)))
    % %        g = ones(size(g)) * 1e6;
    % %    end
    
end





