function [data2, info] = SI2con_2DFLASH_whole(data1, data2, info, maska, skip, precontrast, doT1map, doDynSeq)
% conversion to a value linear w.r.t. concentration
% input args:
%   data1, data2, info - input data
%   skip - how many frames to skip
% output args
%   updated info and data2

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

seqs = length(data1);   % number of precontrast seqs.
slices = info.sections;
info.acq.skip_frames = skip;

fas = 1:90; % flip angle axis

%frames = length(data1{info.acq.kM0Range(1)});
X = size(data1{info.acq.kM0Range(1)}{1},2);
Y = size(data1{info.acq.kM0Range(1)}{1},1);

si = zeros(1,50);
% figure;
% subplot(121)
% hold off;
% plot(fas(1),si(1));
% hold on;


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.FA(seq)*ones(1,frn)/180*pi) ];
    numframes = numframes + frn;
end

img = zeros(Y,X,numframes);   % sequence , frame
ind = 1;
TR = info.acq.TR(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
    indtable(seq,1)=ind;
    for i=1:frn
        img(:,:,ind) = double( data1{seq}{i+skip} );
        ind = ind + 1;
    end
    indtable(seq,2)=ind-1;
    if info.acq.TR(seq) ~= TR
        disp('Error: different TR 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);
    M0m_lowangle_approx = zeros(Y,X);
    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)';
                
                [bla ind_min] = min(info.acq.FA); %... minumum flip angle
                M0_odhad = double(mean(ydata(indtable(ind_min,1):indtable(ind_min,2) ) ) / sin(bla/180*pi));
                
                [E10 M0] = SI2con_fit(xdata, ydata, TR, M0_odhad);
                E10m(y,x)=E10;
                M0m(y,x)=M0;
                
                T10 = -TR/log(E10);
                T10m(y,x)=T10;
                
                M0m_lowangle_approx(y,x) = M0_odhad;
                
                yy = M0 * sin(fas/180*pi) * (1-E10) ./ (1-E10*cos(fas/180*pi));
                figure(10); plot(xdata*180/pi, ydata, 'o', fas,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;
    info.maps.M0_lowangle_approx = M0m_lowangle_approx;
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
    
%%%info.maps.M0  =     M0m_lowangle_approx;

    % if skipmapfitting
    %     load(sprintf(T10mapName, input));
    %     load(sprintf(M0mapName, input));
    %     info.maps.T10 = T10m;
    %     info.maps.M0  = M0m;
    % end
    
    
    flip_angle_dyn = info.acq.FA(end);
    % % get zero signal level:
    % idyn = find(info.acq.FA(info.acq.kM0Range(1):info.acq.kM0Range(2))==flip_angle_dyn);
    % zl = mean(img(:,:,idyn,info.acq.skip_frames+1:end),4);
    
    %load(sprintf(dynamic_seq_mat, slice));
    dyn_frames = length(data2);
    for frame = 1:dyn_frames
        im = data2{frame};
        figure(60);imagesc(im);colormap(gray);
        %                 a = (im-zl)./M0m/sin(flip_angle_dyn/180*pi);
        %                 b = (1-E10m) ./ (1-E10m*cos(flip_angle_dyn/180*pi));
        %                 argl = (1-(a+b)) ./ (1-(a+b)*cos(flip_angle_dyn/180*pi));
        %                 argl(argl<0.00001)=0.00001;
        %                 argl(argl>1)=1;
        %                 A{frame} = -1/info.acq.TR * log( argl );
        nom =    im - info.maps.M0*sin(flip_angle_dyn/180*pi);
        denom =  im*cos(flip_angle_dyn/180*pi) - info.maps.M0*sin(flip_angle_dyn/180*pi);
        argl = nom./denom;
        argl(argl<0.00001)=0.00001;
        im = single(-1/info.acq.TR(end) * log( argl ) - 1./info.maps.T10);
        %%%im = single(-1/info.acq.TR(end) * log( argl ) );
        im=real(im);
        %im(im<0)=0;
        im(~isfinite(im))=0;
        im(~maska) = 0
        data2{frame} = im;
        
        figure(70);imagesc(medfilt2(data2{frame},[5 5]));colormap(gray);title(sprintf('frame: %d',frame))
    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

end
